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Abstract 


la  our  thenaodynaaic  approach^  crack  densities  have  been  treated  as  a 
hierarchy  of  order  paraiaeters  operating  at  successively  srjaller  scales  (froa 
Kscrosccpic  down  to  laicrocrack  scales)  -  The  fornsalisa  renorsializes  surface 
energy  densities  due  to  interaction  between  inicrocracks,  yields  by 
thenaodynamic  self-consistency  effective  elastic  constants  (similar  to  those 
in  effective  medium  approximations)  and  stress  intensity  factors  (SIF)  in 
cracked  solids  and  provides  thermodynamic  definitions  for  observables  (like 
R-curves  and  SIF)  . 

The  fracture  behavior  under  fixed  grip  or  constant  stress  conditions 
has  been  studied  numerically,  above  a  critical  strain  the  material  attains 
an  equilibrium  non-zero  crack^ensity  that  approaches  the  percolation 
threshold  asymptotically  with  increasing  strain.  The  effects  of  cracW^radc 
interaction^  of  anisotropic  crack-density  and  of  pre-existing  pores  have 
also  been  studied.  Crack  propagation  in  ceramic  materials,  subject  to 
microcrack  formation  and  coalescence,  has  been  examined  with  models  that 
differ  in  the  relative  time  scale  involved  in  the  internal  relaxation 
processes  and  the  external  loading  rate.  It  was  found  that  when  the 
external  loading  rate  is  high,  a  rising  R-curve  is  obtained,  which  is 
primarily  due  to  microcrack  shielding.— - — _ 

In  a  numerical  calculation  of  the  Gruneisen  parameters  in  a  cracked 
solid  it  was  found  that  these  decrease  with  increasing  crack  density,  with 
the  rate  of  decrease  dependent  on  the  type  of  effective  field  approximation 
(self-consistent  or  differential  scheme)  that  is  used.  A  Voronoi-polygon 
model  of  granular,  disordered  solids  shows  that  at  fracture  the  number  of 
intergranular  cracks  varies  with  the  line-  or  sample  size  L  as 
Experimental  justification  for  the  critical,  percolation  model  was  found  in 
the  scale  independence  of  the  fragment  size  distributions  in  explosively 
shattered  reinforced  concrete  slabs. 
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PART  1. 


ADHJHISTBAXIVZ 

gtatapwot  of  WorK 

The  purpose  of  this  research  is  the  application  of  percolation  theory, 
information  theory,  molecular  dynamics  simulation  and  fractal  analysis  to  an 
understanding  of  the  basic  issues  in  damage  grovfth,  damage  dependent 
non-linear  mechanical  response  and  fragmentation  processes  in  concrete  and 
rock  under  transient  high-stress  loading.  Although  the  starting  point  will 
be  statistical,  our  ultimate  goal  is  the  merging  of  the  statistical  and  the 
material-oriented  dynamical  approaches  into  a  single  theory.  The  defect 
statistics,  hitherto  regarded  as  homogeneous  in  time  and  space,  will  be 
biased  by  physical  considerations  so  as  to  include  either  process  dependent 
properties  or  specific  material  descriptions.  The  outcome  of  our  theory 
will  be  erstwhile  exemplified  in  a  few  basic,  idealized  cases,  such  as 
isotropic  or  anisotropic  solids,  materials  with  structure  (e.g.,  concrete 
containing  randomly  oriented  ingredients) .  We  shall  investigate  model- 
independent  properties,  including  the  existence  of  critical  points,  the 
nature  of  singular  behaviour,  critical  exponents,  etc.  The  mechanical 
properties  of  cracked  material  near  the  fragmentation  critical  point  will  be 
analyzed  in  terms  of  critical  exponents. 

After  having  laid  down  a  formalism  for  the  micro-crack  based  theory  of 
fracture  in  the  previous  year  of  research  (summarized  in  our  first  Annual 
Technical  Report) ,  we  have  pressed  forward  in  the  original  direction  as  well 
as  in  several  new  directions.  Common  to  all  our  efforts  is  the  statistical 
physios  description  of  fracture  processes  based  on  a  microdefect 
population. 

Since  most  works  in  fracture  are  formulated  in  terms  of  one  or  other 
length  scale,  it  was  important  to  establish  a  consistent  linkage  between 
different  hierarchical  length  scales  (from  the  sample  size  down  to  the 
smallest  relevant  scale) ,  This  was  done,  along  with  other  matters,  in  our 
work  set  out  in  Part  2.  This  formal  work  was  supplemented  by  a  quantitative 
study  (Part  3)  that  showed  how  crack-clusters  affected  mechanical  strength 
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and  fracture  under  sinrolified  but  precisely  defined  conditions.  In  essence, 
we  have  explored  the  elastic  constants  (like  Young' s  nxjdulus)  near  the 
critical  percolation  region-  Curves  of  crack  resistance  (R-curves)  or 
stress  intensities  Kj  have  been  calculated  in  Part  4  to  link  up  the  working 
zone  concept  with  the  halo  of  microcracks  that  form  around  the  propagating 
macrocrack.  The  rising  R-curves,  that  are  found  in  our  calculations  and 
that  have  been  a  puzzle  for  some  time,  are  due  to  the  energy  investment  in 
the  microcracks. 

The  statistical  physics  side  of  fracture  has  its  counterpart  in  the 
probabilistic  treatments  of  fragmentation.  These  are  reviewed  and  developed 
in  Part  5  from  an  information  theoretical  angle.  Worthy  of  emphasis  is 
section  5  where  results  of  percolation  theory  are  incorporated  in  a  Maximum 
Entropy  framework. 

Elastic  moduli  and  Grueneisen  parameters  are  affected  by  the  crack 
population  in  a  stressed  solid.  Their  calculation  requires  averaging  to  be 
performed.  This  is  done  in  Part  6,  which  highlights  the  difference  between 
two  averaging  schemes  (the  self-consistent  and  differential) . 

In  the  last  couple  of  years  the  statistical,  microcrack-based  fracture 
concept  has  gained  strong  following,  though  not  necessarily  in  the  specific 
form  that  is  used  by  the  Soreq  group.  In  a  recent  book  some  of  the 
alternative  approaches  (we  re-emphasize,  each  based  on  a  multiplicity  of 
microcracks  and  employing,  as  a  rule,  Monte  Carlo  methods  and  scaling  laws) 
have  been  sximmarized.  A  review  of  this  book  is  given  in  Part  7. 

Part  8  is  a  preliminary  report  on  how  fracture  develops  in  a  random 
granular  materials.  It  has  a  finite  size  scaling  result  and  shows 
geometrical  correlation  in  the  pictures  of  the  breakage  pattern. 

Some  experimental  vindications  of  the  percolating  multicrack  concept 
appear  in  Part  9,  The  support  cited  is  partly  quantitative,  statistical; 
namely,  through  showing  that  in  different  sets  of  experiments  on  concrete 
the  fragment  statistics  agrees  with  that  of  a  percolating  structure  near  its 
criticality.  Finally,  microphotographs  show  macrocrack  formation  (during 
creep  in  some  ceramics)  as  a  percolative  process  involving  microvoids. 

The  domain  of  applicability  of  our  statistical  physics  picture  for 
rupture  is  not  yet  delineated.  In  particular,  it  is  not  clear  how  it 
depends  on  structure,  composition  randomnsss,  treatment,  the  speed  of 
fracturing  process  or  of  other  parameters  of  the  process.  Theory  and 
cleverly  designed  experiments  will  no  doubt  answer  these  questions  in  the 
near  future. 
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Details  of  Worlc  gerfonned 

Having  found  earlier  an  analytic  law  describing  the  decrease  of 
fracture  strength  with  porosity  and  its  increase  with  the  size  of  pores,  we 
have  shown  that  e.g.  critical  strains  can  depend  markedly  on  such  delicately 
tuned  parameters  as  specific  crack-surface  energy  and  the  anisotropy  in  the 
crack-crack  interaction.  The  history  of  local  stress  near  a  macrocrack 
(namely,  whether  it  relaxes  faster  or  slower  than  the  applied  stress)  was 
found  to  alter  the  computed  R-curves.  The  statistical  spread  in  granular 
structures  (modeled  in  our  work  by  Voronoi  polygons)  makes  the  initial 
breakage  pattern  more  diffusely  located  than  in  a  uniformly  patterned 
granular  system  (as,  e.g.  a  chessboard)  and  makes  percolation  theory  a  more 
realistic  approximation.  The  dynamics  of  flow  in  a  percolating  or  fractal 
framework  has  gained  better  understanding  in  the  last  few  years  and  has  not 
been  investigated  by  ourselves. 

Several  criteria  have  been  used  in  our  research  to  describe  the  opening 
of  a  microcrack  or  the  propagation  of  a  maorocrack,  leading  to  fracture. 
The  Voronoi-method  has  employed  a  fixed  relative  elongation  to  open  an 
interfacial  crack.  This  would  appear  to  favor  the  opening  of  cracks 
separating  small  neighboring  grains;  yet  our  computed  results  show  only  a 
slight  preference  for  this  (Part  8,  Fig.  3) .  The  strongly  rising  R-curve  is 
a  sign  for  the  great  energy  penalty  involved  in  opening  new  microcracks 
around  a  propagating  macrocrack.  Our  computed  results  (Part  4)  show  that 
whenever  this  is  encountered  in  practice,  a  strongly  damaged  working  zone  is 
indicated. 

The  interrelated  aspects  of  micro-crack  initiation,  propagation  and 
coalescence  at  various  hierarchical  levels  of  length  scales  have  been 
treated  in  Part  2.  It  was  shown  how  the  mechanical  parameters  (like  crack 
surface  energy  densities  or  stress  intensity  factors)  differ  at  various 
levels  and  what  is  the  theoretical  framework  (not  yet  worked  out,  at  this 
stage)  that  links  them.  The  theory  also  provides  relationships  betweeri  the 
reduction  in  the  elastic  strengths  (e.g,  in  Young's  modulus),  that 
traditionally  define  "Damage",  and  geometrical  properties  of  crack  clusters 
(namely:  the  effective,  excluded  volume  or  the  cube  of  end-to-end  distances 
(Fig.  1  in  Part  3)].  Plots  of  Young's  moduli  have  been  given  as  function  of 
the  crack  occupation  probability  p  (Part  3,  Fig,  2);  we  have  demonstrated 
that  near  the  critical  defect-  concentration  different  ways  of  averaging 
(i.e.,  different  effective  field  theories)  give  reductions  in  that  are 
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markedly  different.  Some  experimental  efforts  ought  to  be  devoted  to  sort 
out  the  truth. 

It  is  well  known  that  a  fractured  or  fragmented  solid  is  not 
necessarily  at  its  absolute  thermodynamic  energy  minimum,  whence  the 
difficulties  in  the  application  of  thermodynamics  to  fracture.  On  the  other 
hand.  Maximum  Entropy  Methods  (MEM)  are  based  on  states  of  knowledge  (which 
are  at  the  roots  of  all  probabilistic  reasoning) ,  which  also  hold  good  at 
non-equilibrium.  In  Part  5  we  have  given  MEM  predictions  for  fragment 
statistics  both  for  a  high-strain  rate  single,  sudden  fracture  process  and 
for  (still  a  high  strain  rate)  situation  in  which  fracture  develops  over  a 
spread  of  time.  Though  the  statistics  is  for  fragments,  it  is  the 
proliferation  of  cracks  that  forms  fragments. 

Predictions  of  our  statistical  physics,  micro-crack  approach  depend  on 
the  values  of  a  relatively  large  number  of  parameters  and  on  the  physical 
conditions  (experimental  situations)  under  which  theory  is  tested.  Therefore 
a  fully  convincing  validation  of  our  approach  appears  possible  only  through 
the  performance  of  a  set  of  planned  experiments  under  controlled  conditions 
(and  not  the  interpretation  of  published  experimental  results) .  Still,  two 
important  facts  of  our  tneories  have  received  direct  experimental  support 
(Part  9)  :  (1)  the  merging  into  macrocracks  of  smaller  defect  units 
(microvoids)  has  been  documented  under  slow  degradation,  high  temperature 
conditions  (creep) ,  (2)  quantitative  evidence  has  been  found  from  medium  and 
large  scale  deformation  experiments  on  (reinforced)  concrete  for  our 
prediction  (in  1985)  that  fragmentation  is  a  critical  process  due  to  the 
approach  of  the  crack  density  to  a  critical  value  (Pen)  •  observed 
mass-distributions  of  fragments  in  two  sets  of  experiments  follow  quite 
closely  the  percolation  theoretical  predictions. 

Finite  size  scaling  (i.e.  dependence  on  sample  size)  of  the  number  of 
broken  springs  has  been  derived  in  Part  9  using  a  computer  code  that  first 
generates  a  random  structure  and  then  follows  its  time  development  under 
applied  stress  or  strain. 
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la  tisa  tvo  years  of  research  tiie  rsicrccraci:  pi»3latioa  theory  was 
ceyelcred  to  interpret  fractcre.  The  basic  theory  was  prodaced,  cateriai 
rcptare  behavior  was  dsscriS^d  by  calcalatiag  crack  resistance  -  (R)  curves 
aaa  ezperireatal  ccafiraatica  was  established  froa  fragaeat  statistics  in 
ccacrete  and  fsca  rsicro-phctographs  of  void  coalescence  in  the  creep  of  sors 
ceraaics.  «e  have  derived  frca  oar  crack-cluster  description  a  sharply 
rising  R-carve  frequently  seen  in  ezperiaent.  Observed  statistics  cirisiag 
frcn  fracture  (as  in  fragreatatioa)  obey  a  power  law,  in  accordance  with  the 
critical  behavior  occurring  near  the  percolation  limit. 


9 


psar  2 


caacs-CLUSTsa  dbscsiptios  o?  paacsnas 


T  sHERHODYKsmcai.  sssaEf  ?micriCNai.s 


R.  Snglman,  M.  Murat  and  Z.  Jaeger 
Soreq  Nuclear  Research  Center 
Yavne  70600,  Israel 


Abstract 

The  traditional  physics  approach,  as  exemplified  by  lattice  dynamics, 
can  be  used  to  develop  a  thermodynamic  formalism  that  generalizes  continuum 
mechanics  (elasticity  theory)  by  including  a  hierarchy  of  cracks,  namely 
microcracks  that  join  up  to  form  macroscopic  cracks .  Fracture  processes  are 
derived  by  minimizing  thermodynamic  energies  that  are  functionals  of  crack 
densities  and  of  strains  (or  stresses)  that  vary  on  successively  smaller 
scales  (from  macroscopic  down  to  microcrack  scales) .  The  formalism 
renormalizes  surface  energy  densities  due  to  interaction  between 
microcracks,  yields  by  thermodynamic  self-consistency  effective  elastic 
constants  (similar  to  those  from  effective  medium  approximations)  and  stress 
intensity  factors  (SIF)  in  cracked  solids  and  provides  thermodynamic 
definitions  for  observables  in  engineering  practice  (like  R-curves  and  SIF)  . 
Crack  densities  can  be  treated  in  the  formalism  as  order  parameters: 
criteria  for  fracture  using  this  approach  are  developed  in  the  sequel 


paper. 


-  2  - 
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1 -  Introduction 

A  theraodynamic  theory  of  lattice  dynamics,  in  which  interatomic  bonds 
are  the  discrete  units,  was  con5>lete  by  the  1930^s  ^  .  The  continuum  theory 
of  elasticity  not  based  on  atomistic  units  dates  back  to  the  days  of 
Rayleigh^' To  describe  fracture  a  need  has  arisen  for  an  intermediate 
type  of  theory  of  solid  behavior  in  which  there  are  both  discrete  entities, 
namely  cracks  of  various  sizes,  alongside  with  intact  pieces  of  the  material 
which  can  still  be  represented  by  continuum  mechanics.  There  is  a  good  c  'al 
of  experimental  support  for  using  this  approach  to  explain  fracture  in 
either  brittle  or  ductile  solids.  The  argument  is  stated  in  more  detail  in 
^3.5  of  this  work  with  some  basic  referencing  and  is  intended  to  justify 

C 

the  crack-cluster  origin  of  fracture  in  brittle  materials  and  the  microvoid 

c  n 

coalescence  theory  of  ductile  failure  '  . 

When  looked  at  from  a  broad  historical  viewpoint,  this  paper  and  its 

O 

sequel  can  be  regarded  as  an  example  of  the  recent  trend  to  give 
engineering  practices  in  fracture  mechanics  a  precise,  physics  oriented 
formulation.  Characteristically,  both  thermodynamics  and  some  form  of 
statistical  mechanics  are  brought  into  play:  the  first  in  the  form  of 
employing  internal  variables  and  their  conjugate  forces  and  the  second  by 
relating  macroscopic  properties  of  the  material  to  microscopic  variables. 
Both  features  were  already  present  in  Rice's  treatment  of  metal  plasticity 
in  terms  of  slip  kinetics^' . 

The  more  recent  efforts  that  form  part  of  the  trend  are  too  numerous  to 
list  exhaustively.  One  wishes  to  mention  Krajcinovic  et  al.'s 

derivation  '  of  the  conjugate  force  to  the  Budiansky  -  O'Connell  damage- 

.  1^ 
volume  variable  . 
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She  present  thermodynamic  app^roach  differs  from  previous  works  in 
several  respects.  First,  our  approach  is  presented  as  a  hierarchical 
extension  of  Born-Huang's  classic  treatise  of  lattice  dynamics,  when  applied 
to  discrete  units  (cracks)  of  increasing  sizas-  Admittedly,  this  is  only  a 
didactic  device  but  it  facilitates  the  systematic  derivation  of  results. 
Next,  the  transition  between  micro  and  macrodefects  is  performed,  smoothly 
by  the  use  of  crack  clusters  of  increasing  size,  in  a  way  originated  by 
percolation  theory.  The  energy  release  term,  also  called  the  mechanical 
energy,  is  defined  for  crack  clusters  in  terms  of  the  stress-free  volume 
that  is  the  equivalent  of  the  Budiansky  -  O'Connell  damage  -  volume  for 
simple  (e.g.  penny-shaped)  crack.  We  then  show  how  several  effective  medium 
approximations  can  be  derived  from  simple  consistency  requirements  of  the 

energy  functionals.  Further,  we  include  crack-crack  interaction  in  the 

energy-functional,  formally  in  a  manner  natural  to  many-particle  solid-state 
physics,  but  not  yet  met  with  in  fracture  mechanics.  We  indicate  how  the 

interaction  leads  to  a  renormalization  of  observables  (like  crack  surface 

energy  densities  or  stress  intensity  factors)  in  a  way  that  can  account  for 
the  large  numerical  discrepancies  between  quantities  that  correspond  to 
micro-  and  macrocracks.  In  addition,  we  propose  a  thermodynamically 
acceptable  definition  of  crack  resistance  -  (R-)  curves  and  of  other 

observables  on  which  the  engineering  characterization  of  materials  rests. 
Finally,  we  introduce  crak  densities  as  order  parameters  in  the  energy  and 
obtain  fracture  mechanics  by  minimization  of  the  energy  with  respect  to 
crack  densities.  This  extension  of  the  (more  generally  valid)  theory  is  the 
basis  of  the  numerical  study  in  the  sequel  paper®.  The  defect  density  also 
appeai:s  as  a  free  variable  in  the  thermodynamic  theory  of  Gelbart  and 
cowork«irs^^ . 
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The  present  approach  is  not  free  of  uncertainties,  which  v;e  try  to 
describe  as  we  go  along.  However,  we  feel  that  these  have  been  in  the 
subject  for  some  time.  Thus,  the  definition  of  R-curves  is  notoriously 
ambiguous^^'^^.  The  validity  of  any  effective  medium  description  has  been 
questioned  for  a  medium  containing  microcracks  .  The  demarcation  between 

,  n 

micro-  and  macrocracks  is  fuzzy  .  Yet,  we  feel  that  the  present  formalism 
can  be  helpful  as  long  as  its  limitations  are  clearly  recognized. 

2.  The  Energy  Functionals 

Our  objective  is  a  theory  that  is  intermediate  between  continuum 
elasticity^' ^  and  discrete  lattice  dynamics^' in  the  sense  that  the 
displacement  variable  e(r)  (r  =  the  position  in  the  material)  is  decomposed 
into  a  macroscopic  strain  and  local  strains  ej^(Arj^)  around  the  i 

discontinuity  (crack,  or  in  ductile  solids,  void)  in  the  damaged  solid.  Ar  = 
r  -  rj^  where  r^  is  the  location  of  the  i'th  crack.  The  distribution  of  the 
cracks  in  the  solid  p(r)  can  be  preassigned  or,  else,  can  be  treated  as  the 
independent,  order-parameter  in  the  energy  functional 

U  =  U({e(r)  },  {p{r) )) 

The  latter  approach  is  described  in  ^3.5  and  is  developed  in  the  sequel®. 
Following  standard  lattice  dynamics  treatments  ',  we  identify  with  U  the 
internal  energy,  which  is  written  as  a  sum  of  integrals  over  the  full  volume 
V  of  the  material,  which  includes  several  cracks.  (So  as  to  put  over  our 
ideas  with  the  best  effect  we  dispense  with  the  tensorial  indices  of  strain 
(or  stress),  which  practice  is  frequently  followed  in  the  schematic 
representations  of  fracture  processes  .  In  a  quantitative  treatment, 
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with  the  precision  of  e.g.  Ref.  13,  our  formalisin  should  take  account  of 
tensorial  properties,  of  non-zero  Poisson  ratios  V,  etc.]  We  write 


1 

U(8)  =  -  Eo 

1 

(r)dr  -  —  EpI 

d(Ar)  +  Yo  ^ 

2 

2  i  . 

i  . 

Ar  e  5^  ^ 


+  I 
ij 


d(Ari)d{Arj)  DCr^-rj) 


(1) 


The  first  term  is  the  usual  strain  energy  of  the  solid.  The  second 

term  gives  the  reduction  of  the  energy  due  to  the  relaxation  from  a  strained 

state  in  a  volume  around  the  i'th  crack.  (  is  a  representative 

ft 

end-to-end  distance  of  the  crack”.)  This  energy  release  is  the  familiar 
"mechanical  energy"  appearing  in  the  fixed  stress  configuration  in  the  form: 


1 

-  -  (2) 

2 

(Ref.  4,  p.  253) .  The  exact  form  of  the  stress-free  volume  depends  on  the 
crack's  shape  and  orientation  and  requires  a  more  careful  treatment  than  can 
be  given  here;  but  there  can  be  no  doubt  about  the  relevance  of  this  term  to 
fracture.  Thus  this  term  is  responsible  in  the  Griffith  theory  for  crack 

instability,  that  comes  about  when  the  macroscopic  stress  (r„  exceeds  the 

1  ft  1  0  . 

(tensile)  fracture  stress  cr^  '  .  The  meanings  of  the  Young  moduli  and 

Eg  in  (1)  and  of  E^^  in  (2)  will  be  presently  explained,  in  the  course  of  the 

derivation  of  the  effective  elastic  constant  Egff  in  the  cracked  medium. 

The  third  term  represents  the  surface  energy  needed  to  open  up  the 

discontinuities  in  the  zeroth  order  (non-interacting  microcracks) 

approximation  with  the  basic  surface  energy  density 
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The  last  term  is  the  crack-crack  interaction.  It  is  decomposed  into 
two  parts;  a  strain  independent  and  a  strain  dependent  part  in  the  form 

D(rjL-rj)  =  (7(r^-rj)  +  Eg  (Ar£)ej  (Ar j)  K(rj^-rj)  (3) 

There  is  no  linear  term  in  e  since  both  the  macroscopic  strain  and  the 

local  strains  e£(Ar£)  are  zero  in  the  equilibrium,  unstressed  configuration. 

[The  argument  is  sinilar  to  that  leading  to  the  elimination  of  the  linear 

terms  in  the  atomic  displacement  from  the  lattice  potential  of  Born  and 

Huang^,  III  ^11) .  Under  conditions  of  fixed  temperature  the  free-energy 

function  F  would  be  used  instead  of  U.  This  differs  from  U  by  a  crack- 

21  14 

configuration  entropy  term  that  has  been  considered  by  us  '  and  others'^’. 
However,  it  is  not  an  essential  component  in  the  theoretical  development  and 
will  be  disregarded  for  the  present.  More  important  is  the  entropy  change 
AS  associated  with  the  irreversible  breaking  of  the  bonds.  This  will  be 
treated  in  §3.4. 

The  energy  function  U  leads  naturally  to  a  formal  treatment  of  a  number 
of  effects  that  are  central  to  cracked  crystals.  We  list  these  now  and 
proceed  to  their  discussion  one-by-one,  emphasizing  that  in  a  complete 
treatment  all  the  effects  enter  together:  we  separate  them  for  reasons  of 
clarity. 

(I)  Self  consistent  effective  elastic  constants  from  thermodynamical 
considereations . 

(II)  The  macrocrack  surface  energy  density  Y{^'Yo)’ 

(III)  Self  consistent  stress  intensity  factors  (SIF)  at  tips  of 

interacting  cracks. 

(IV)  Thermodynamic  definition  of  R-curves,  stress  intensity  factors 

(V)  Crack  densities  p  as  order  parameters. 
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3,1  Effective  elastic  moduli 


To  treat (I)  we  simplify  (1)  as  follows: 

(a)  Make  U[8(r)]  the  function  of  one  single  strain  variable  only,  this 
ignores  all  strains  (Arj^)  associated  with  cracks  (or  cracks  tips)  except 
as  averages. 

(b)  Consider  at  the  first  stage  only  non-interacting  cracks;  so  that  in 
(1)  the  last  term  is  ignored. 


Then  the  strictly  dependent  part  of  U(eeo)  has  the  form 


Eg 

-  Eg  - E^  V  [1 - 

2  2  Eq  V 


yex 


where  V  is  the  volume  of  the  solid  and  V®^  is  the  "excluded"  crack  volume, 
given  by 

i 

analogous  to  Budiansky  and  O'Connell's^^  <a^b>.  We  emphasize  that  V®^  is 
non-zero  even  if  the  cracks  are  perfectly  flat  and  thin,  i.e.  do  not  occupy 
any  volume  in  the  solid.  This  being  so,  we  can  transform  the  constant  grip 
internal  energy  in  (4)  to  a  constant-stress  energy-functional  H(cf„) 
analogous  to  the  enthalpy®,  by  a  simple  Legendre  transformation  not 
involving  V®^  (since  no  work  is  done  by  to  separate  the  flat  cracks)  . 
Therefore 


H  (cr„)  =  u  [e„  (a„)  ]  -  (r„e„v 


(5) 


-  8  - 


16 


Actually  our  starting  point  might  also  have  been 


H(cr„) 


1  ff- 


2  E, 


1  (T„ 


V  -  - 


,ex 


2  E. 


(6) 


To  relate  Eg  in  (4)  to  E^^  in  (6)  we  use  the  thermodynamic  definitions  of 
conjugate  variables,  namely: 


1  0H(cfJ 

- =  -  e„  (7) 

V 


1  Ou((r„) 

-  =  (8) 

V 

and  the  definitions  of  the  effective  Young's  module  in  the  solid 

containing  cracks  (which  for  the  moment  regarded  as  non-interacting) 


U(eJ 


2E, 


eff 


(9) 


(10) 


The  presence  of  the  same  E^gg  in  (9)  and  (10)  amounts  to  the  self- 
consistency  requirement,  emphasized  by  Hashin*''  and  by  other  workers 
therein.  The  consequence  of  self-consistency  can  be  put  in  the  form  (which 
is  not  unique) 


Eeff  =  EqCI  +  (0-l)v]/(l+9v) 


(11) 
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Eg  =  Eg/(l+0v) 

(12) 

Ejy  =  Eg[l-(l-0)v] 

(13) 

where  we  introduce 

V  H  v®^/V 

and  0  is  some  quantity  between  0  and  1.  The  limiting  case 

0  =  0;  Egff  «  Eq(1  -  v);  Eg  =  E^;  E^^  =  (14) 

.  .  1  ^ 

13  the  choice  adopted  by  Budiansky  and  O'Connell  /  characterized  by  an 
effective  Young's  module  that  goes  to  zero  at  v  =  1,  equivalent  to  ref 
(13) 's  9/16  catastrophy.  At  the  opposite  extreme 

0=1;  Egff  =  Eg/(1  +  V);  Eg  =  Egff;  E^  =  E^  (15) 

®eff  tends  to  zero  only  asymptotically  for  high  crack  density  as  also  in  the 
Differential  Scheme  (DS)  of  Hashin^^. 

It  is  important  to  recognize  that  the  class  of  effective  moduli  in  (11) 
-  (13)  were  obtained  by  thermodynamic  self-consistency  and  not  by  any 

physical^^  or  geometrical^^  argument. 

Extension  of  these  considerations  to  account  for  the  tensorial 
character  of  E  and  cx,  including  effective  Poisson's  ratios  Vg£f  will  be 
given  in  a  future  publication. 

3.2  Renormalized  Surface  Energy  Density 

The  numerical  value  of  the  surface  energy  density  y  varies  according  to 
whether  the  created  discontinuity  is  on  an  atomic  scale/  micro  or  macrocrack 
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scale  ' '  .  The  breaking  of  atomic  bonds  is  too  energy  costly  (~  10  eV) 
to  be  caused  by  ordinary  fracture  processes;  7  -  values  for  visible 
discontinuities  are  lower  than  these  by  a  factor  of  10~^  -  10”^  [Ref.  26]; 
the  presence  of  microcrack-clouds  immediately  surrounding  the  macrocracks 
and  creating  the  "plastic  zone"  near  the  crack  tip  reduces  ^^9^27,28^ 

The  first  term  in  (3)  is  the  source  of  these  effects.  It  represents 
the  mutual  influence  of  the  i'th  and  j'th  cracks.  Accordingly,  the  Ar 
integrations  in  the  last  term  in  (1)  are  to  be  taken  over  the  discontinuity 
surfaces.  Regarding  as  atomic  quantities  of  the  order  of  dissociation 
energies  the  integral  in  (1)  is  to  be  taken  over  all  cut  atomic  bonds  in  the 
microcracks.  For  the  i'th  microcrack  the  last  integral  is  then  written  as  a 
correction  to  the  atomic  quantity  y^  in  the  form: 


SYi 

1 

dSi  -  ti 

' 

dSi  1 

d(Arj  ) 

V(ri  -rj 

) 

(16) 

1 

1  1 

1 

jl=il 

1 

1 

1 

1  the 

integrations 

dS,d(Arj  ) 

are  over  atomic 

dimensions. 

tj  is  the 
*1 

(atomic)  thickness  of  the  ii  microcrack  and  r.:  is  the  coordinate  of  the 
jl(=il)'th  atomic  thickness  of  the  i]^  microcrack.  In  brief,  and  in  a 
schematic  form,  as  given  by  (16)  represents  the  correction  to  y^  due  to 
other  broken  atomic  bonds  in  the  i^  microcrack: 


Yi  =  Yq  +  5Yi 
1 


1 
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For  a  crack  cluster,  consisting  of  a  connected  configuration  of 
several  microcracks,  and  denoted  by  i2/  an  integral  similar  to  (16)  will 
further  renormalize  ,  as  follows: 


7i  =  Yi  +  SYi 

2  12 


where  Sy^  represents  contributions  to  the  i^  microcrack  by  other 

microcracks  belonging  to  the  i2  crack  cluster. 

Microcracks  from  other  crack  clusters  unlikely  to 

contribute  to  Yi  since  the  interaction  O'  is  of  short  range.  However, 


close-by  microcracks  surrounding  the  i2  crack  cluster  will  alter  Yi  /  I*'  the 

•  19 

sense  of  the  plastic  zone  correction  to  Y  due  to  Irwin 

While  the  preceding  discussion  is  overly  qualitative,  it  can  lead  the 
way  to  a  systematic  prediction  of  surface  energies  observed  in  fracture. 


3 . 3  Self  consistent  stress  intensity  factor  (SIF)  at  crack  tips . 


We  now  go  beyond  regarding  the  energy  functional  U  as  depending  on  only 
one  variable,  the  macroscopic  strain  only.  Various  extensions  are 
possible .  Indeed  we  shall  not  go  along  the  whole  length  with  Born  and 
Huang^  to  de..oribe  U(E)  as  a  function  of  all  atomic  displacements.  Rather, 
we  consider  strains  only  in  the  vicinity  of  cracks.  Here,  too,  we  can  have 
strains  e{Ax:^)  associated  with  the  i  crack  as  function  of  all  points 
near  the  crack,  or  only  on  tne  crack  surface.  Leaving  the  latter 
possibility  aside,  we  shall  describe  more  simply  (and  less  accurately)  local 
crack-induced  strain  effects  in  terms  of  a  single  strain  associated  with 
each  crack  (or  crack  tip)  .  We  need  not  make  at  this  stage  a  distinction 


betwses  sicrocra'Cis,  crzzk.  clcsters  cr  assrscrssSa.  ISag  ffsiarnlstasa  wSli 


be  sboTO  to  be  es^valeint  t©  srerfass  c?s'^C3?st£fflr:s  cf  tbe  effectice  Sir  at 
crack-tips  &;2  to  iEteractirg  cracks  ©r  at  a  aacrocrack  dee  to  sererai 
rdcrocracks,  or  a  eiacrocrask  array.  Sbe  lerei  sf  la 

Kachanov's  vorks^^'**^  is  sfssilar  to  the  cire^nt  ore.  Stare  is  ks  ditticalty 
in  refining  the  present  sirplified  theory,  czevidsd  it  is  accosted  that 
deeper  levels  of  apprezirratiens  can  be  get  sinply  by  sraking  Sjej  a  fnretion 
cf  additional  strain  variables-  Additional  stress  variables  are  obtained 
the  conjugate  relations. 

1  3o(e) 

- =  ill) 

V 

1  3H((r) 

- -  e.  (18) 

V  da^ 

Inserting  into  the  first  two  integrals  in  (1)  for  regions  near  the  i 

crack 

8(r)  =  £„  + 

where  =0, 
i 

we  obtain  integrands  in  the  form 

+  2e„ei  +  (19) 

where  the  first  term  gives  the  macroscopic  strain  energy,  the  second  term  is 
an  interaction  whose  physical  meaning  is  the  strain  near  the  crack  induced 
by  the  macroscopic  strain  and  the  third  term  is  a  local  strain  energy.  In 
the  crack  interaction  term  the  strain  dependent  part  [shown  in  (3)3  has  the 
similar  forms,  multiplying  the  interaction  K^j, 
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120} 


la  tbs  siseace  cf  sracrasscaiic  strain  tire  SI?  can  be  calcolated  frsra  (17), 
«aisb  giv?es 


^  £3*3 


(21) 


-  ^«"%3  "  %3>3£j 


(22) 


-,3  -3 


(23) 


having  used  the  srraaticn  ccaveatica  for  repeated  cracb-iadices  in  the  last 

<j|> 

two  eaxjressicas  aad  iatrodrciag  the  effective  rxdal'us  aatrix  .  Asscaing 
the  existence  of  the  inverse  *  we  rewrite  the  interaction  in  (1),  (3) 

as 


J  C •  s!d  rC*  t  *  t  ^  “ 

1  I'l'  "  ”3 

i#i'/3/3  11'  3' 3 


(24) 


The  nvean  stresses  on  each  crack  i  being  given  as  in  (23)  we  can 
calculate  the  stress  fields  due  to  them  at  each  point  r^^  on  the  i^^th 
crack,  following  standard  methods  in  elasticity  .  From  these  then  follow 
the  SIF- 

OQ  Oft 

The  result  is  equivalent  to  Kachanov's  '  (after  making  the  necessary 
adjustments  in  the  formalism,  like  V®^  =0,  Eg  =  E^) .  It  is  significant 
that  the  self  consistency  of  the  mean  stresses,  obtained  by  him  in  the  form 
of  a  mean  field  assunption,  again  arises  from  the  thermodyntmic  formulation 


(Eq.  (17),  (18)]  automatically.  Unfortunately,  the  advance  in  the  theory  is 
but  formal  since  crack-crack  interaction  stresses  have  been  calculated  only 
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fcr  a  few,  si^ciai  gecretries^^'^^-  (Sfee  iaterpolatioa  fozzsila  used  in  the 
sequel  paper  has  an  illustrative  purpose.) 

3.4  gheraodynaain  observables 

5he  nacroscopic  variables  0^.,  defined  in  5qs.  (7),  (8),  are  evident 
observables.  So  are  the  local  variables  e-  (Eqs.  (17), {18)3#  though  a^, 
ray  not  be  externally  controllable.  However,  our  free  energy  contains 
additional  pararaters  which  should  also  give  rise  to  observables,  nainely  the 
crack  positions  or  crack  densities.  In  principle  each  microcrack  should  be 
associated  with  a  (cicroscopic)  obseirvable.  From  a  practical  viewpoint  we 
consider  a  much  smaller  number  of  variables,  which  will  be  further  reduced 
to  two:  the  critical  stress  intensity  factor  Kj^;.  and  the  R-curve  variable. 

To  obtain  Kjq  from  the  energy  functionals  U  or  H[{1)  or  (5)]  we 
associate  with  the  crack  cluster  voliur.es  (belonging  to  the  i'^th 

macrocrack)  a  length  a^. 


while  for  a  penny-shaped  or  wedge-like  macrocrack  other  relationships  will 

20  21 

hold.  We  have  already  discussed  at  some  length  in  earlier  works  '  that 
the  dependence  of  on  the  number  of  microcracks  belonging  to  the 

crack-cluster  is  related  to  the  fractal  dimension  of  the  macrocrack  and  only 
in  very  special  cases,  like  extremely  correlated  crack  growth,  will  depend 
on  the  extent  of  the  macrocrack. 

With  certain  assumptions  about  the  general  mechanism  of  crack  growth, 

1ft 

which  are  basically  similar  to  those  of  Griffith-^®  and  will  be  elaborated  on 
in  the  next  section,  we  can  regard  the  a^  as  free  parameters  in  the  energy 
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functional  and  obtain  their  value  a^  at  equilibritm  for  a  fixed  tensile 

T 

stress  d  from 


3H{(y) 

- =  0  (j=all  loacrocrack  labels)  (26) 

while  for  fixed  strain  is  found  from 

Du(ej 

-  =  0  (27) 

3.3 

The  form  of  U  in  (1)  makes  it  plausible  that  for  a  sufficiently  high 
stress  an  accessible  minimum  (local  or  absolute)  value  of  will  arise 
which  will  be  of  the  order  of  the  sample  size:  the  sample  will  then  be 
fractured.  The  least  stress  achieving  this  condition  is  the  fracture 
stress, 


c 


T 


(28) 


The  phrase  "accessible  minimum  (local  or  absolute) "  will  be  amplified 
in  the  next  section;  "minimum"  refers  to  the  energy  as  function  of  the 
microcrack  density  .  The  minimum  changes  depending  on  what  conditions, 
e.g.  fixed  stress  or  strain,  hold,  the  accessibility  of  that  minimum  hinges 
on  the  issue  whether  microcrack  opening  and  healing  occur.  The  critical  SIF 
will  then  be  given  by 


^IC 


(29) 


-  16  - 


24 


where  b  is  the  linear  dimension  of  the  microcrack  and  a  numerical  factor. 
(In  a  percolation  picture  p^.  is  the  critical  crack  occupation  at 
percolation . ) 

R-curve  measurements  are  aimed  at  determining  the  energy  input  rate 
needed  to  extend  a  stable  macrocrack  (not  yet  cutting  across  the  sample)  by 
unit  length.  Techniques  of  measurements  have  been  described  in  detail^^”^® 
and  involve  first  the  application  of  stress  leading  to  a  small  extension  of 
a  crack  and  subsequently  the  measurement  of  the  change  of  the  Young's  module 
(or,  of  its  inverse,  the  compliance)  due  to  the  crack  extension.  Important 
for  our  considerations  is  the  circumstance  that  the  crack-extension  is  an 
irreversible  process,  involving  a  finite  entropy  change  AS,  whereas  the 
Young's  module  is  measured  (before  and  after  the  crack  extension) 

OO  , 

isoentropically  AS  =  0.  As  described  lucidly  by  Wannier  ,  in  hard 
materials  the  irreversible  change  takes  place  almost  without  any  change  of 
stress,  it  is  therefore  permissible  to  write  the  energy  functional  of  (1) 
and  (4)  as 


H((r) 


1 

2  j 


dAj)a 


(30) 


u(e)  “ 


1 

-  (e^/EQff)g 


2 


+  (y^ 


^j>e 


(31) 


where  the  subscripts  are  the  commonly  used  specifications  of  the 
thermodynamic  quantities  which  are  held  constant  in  the  term,  y  is  the 
renormalized  surface  energy  density  in  accordance  with  our  discussion  in 
^3.2,  and  the  integrals  and  sums  Zj  the  total  area  of  all 

microcracks.  The  dependences  on  the  macrocrack  sizes  aj^  enter  in  the  first 


term  through  Eg^f  and  in  the  second  term  through  the  total  crack  area  which 
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is  evidently  a  functional  of  the  linear  macrocrack  sizes  a^^.  The 
irreversible  energy  change  has  been  identified  above  as  this  second  term  and 
its  derivative  as  the  crack  resistance  for  the  i  microcrack 


ba^  j 


(4  =  0,  e) 


However,  in  equilibrium  for  each  macro  crack  size,  when 


dH 


and  similarly  for  U  (e)  . 

It  follows  therefore  from  (32)  and  (33)  that 


(33) 


R 


=  (Eeff)“^  (34) 

2  ba 

=  G(the  stress  energy  release  rate) 

where  the  subscript  i  has  been  dropped  since  measurements  are  usually  made 
upon  a  unique  (usually,  prenotched)  crack.  An  equivalent  relation  holds 
under  fixed  grip  (e  constant)  conditions. 


The  staleness  of  the  result  (34)  is  misleading:  In  actuality,  it 
establishes  R  as  a  proper  thermodynamical  observable  whereas  in  all  previous 
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approaches  known  to  us  R  features  as  the  a  derivative  of  part  of  the  energy, 
which  is  not  a  thermodynamical  observable. 

3.5  Crack  densities  p  as  stochastic  variables 

In  writing  out  the  energy  functional  U  (Eq.  (1))  the  origin  of  the 
cracks  was  left  unspecified  and  so  was  the  precise  nature  of  and 
distinctions  between  macrocracks  (or  crack-clusters),  microcracks  and  broken 
interatomic  bonds.  These  issues  will  now  be  addressed,  against  the 
background  of  considerable  uncertainties. 

The  need  for  a  discussion  arises  because  of  the  growing  realization 
that  discontinuities  in  solids  visible  to  the  naked  eye,  ranging  from 
fracture  to  spall  planes  to  pores,  arise  from  the  merging  or  coalescence 
of  microcracks  (not  of  visible  size)  '  .  Moreover,  there  is  clear 
evidence  for  the  role  of  microcracks  that  surround  a  fracture  crack  in 
toughening  and  in  the  enhancement  of  the  surface  energy  through  the  presence 
of  a  plastic  zone  (we  refer  to  our  discussion  of  the  renormalized  7  and  to 
References  19,27,37). 

In  their  discussion  of  crack  initiation  Curran  et  al.  define 
microcracks  as  those  having  sizes  comparable  to  "graininess"  that  defines 
the  continuum  limit  of  the  material  and  supply  a  list  of  them  to  which  we 
refer  for  factual  details .  The  definition  assigns  to  microcracks  a  measure 
of  stability,  at  least  as  regards  to  their  tendency  only  to  grow  not  to 
closeup.  Opposing  this  viewpoint  is  the  suggestion  that  physical  processes 
of  a  given  type  (such  as  having  a  given  loading  rate)  result  in  microcracks 
of  a  definite  average  size,  implying  the  possibility  that  crack-sizes 


.  .  .  .  40 

fluctuate  also  in  time  until  converging  upon  an  average  . 


Thus,  in 
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essence,  crack-  growth  arises  only  through  a  proliferation  of  cracks, 
through  merging,  as  already  noted  with  respect  to  macrocracks^® . 

For  clarity,  let  us  remark  that  in  this  work,  and  presumably  in  many 
other  modern  papers,  macrocrack  is  taken  as  synonomous  to  (micro-)  crack 
cluster.  The  reluctance  of  macrocracks  to  disappear  by  closing  up  under 
normal  laboratory  conditions  is  apparent. 

There  might  well  be  ^nb-microcrack  entitities  which  could  represent 
sites  of  nucleation  for  microcracks  and  which  open  and  close  until  they 
reach  microcrack  dimensions.  This  hypothesis  is  physically  sensible  since 
extremely  high  stresses  can  concentrate  upon  tiny  surfaces,  and  these  are 
capable  of  changing  the  atomic  configurations  one  way  or  the  other. 

In  summary,  in  addition  to  the  picture  of  a  given  microcrack  population 
which  can  increase  by  further  application  of  stress,  there  are  other 
possibilities  of  fluctuating  populations  (of  microcracks  or  of 
sub-microcracks)  which  will  adjust  to  the  global  criteria  of  thermodynamics. 
In  other  words,  one  may  venture  to  hypothesize  ergodicity  among  the 
different  realizations  of  crack  populations  in  a  solid  that  is  (or  is  not) 
subject  to  applied  stresses.  This  situation  is  amenable  to  a  Monte  Carlo 
simulation  of  crack  populations,  e.g.  by  percolation  techniques. 

To  include  a  stochastic  crack  population,  the  needed  concepts  are  a 
density  p(r)  of  potential  crack  sites,  that  is  a  function  of  the  coarse¬ 
grained  position  r  (so  that  p(r)  varies  slowly  over  ranges  of  several 
microcrack  sizes),  and  the  variables  p(r)  (taking  values  between  0  and  1) 
which  gives  the  probability  of  occupation  of  potential  crack  sites.  Though 
not  absolutely  necessary,  p(r)  is  regarded  as  a  free  variable  of  the  energy 
functionals  U  or  H,  such  that  these  are  minimal  for  the  actual  occupation 
probabilities  p(r); 
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-  or  -  =0  at  p  =  ^.  (35) 

8p(r)  8p(r) 

The  structure  of  the  energy  as  function  of  p(r)  and  of  p{r)  is  given  in 
the  following: 


U({e(r) ),  {p(r) ))  =  -  Eq 

2 


+  Yob" 


1 

e^(r)dr  -  -  Eg 

drp(r) 

2 

. 

•  drp(r)(p(r)  + 

cb2  H?p|2] 

d(Ar)p(r)e^(r) 


Ar6b(r) 


+ 


drdr'p(r)p(r'  )p(r)p(r'  )D(r-r') 


(36) 


Here  b  is  the  linear  size  of  the  basic  crack  unit  (plaquette)  and  a  gradient 
term  (with  a  positive  phenomenological  coefficient  c)  has  been  included  in 
the  third  term  to  discourage  spatial  fluctuations  in  p{r). 


4 .  Discussion 


Some  arguable  assumptions  of  the  theory  have  already  been  discussed  in 
the  article.  One  might  also  cavil  at  the  use  of  flat  cracks,  while  in  other 
approaches  central  place  is  given  to  the  crack  width^^  or  to  crack  opening 
displacements  .  It  is  possible  to  replace  in  our  formalism  [e.g.,  Eq.  (36)] 
the  probability  field  p(r)  by  a  width  "distribution  without  much 
modification.  We  also  note  that  the  macro-cracks  or  crack  clusters  possess 
widths  due  to  their  geometrical  structure. 

Among  possible,  experimental  verifications  of  the  theory  one  would  name 
determination  of  the  microcrack  distribution  near  a  macrocrack,  as  a 


function  of  distance  and  orientation.  This  is  technically  feasible  now^^. 
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Independent  measurements  of  E^ff  and  (the  excluded  volume  due  to  crack- 
clusters)  would  test  the  relations  proposed  in  ^3.1.  R-curve  measurements 
under  controlled  conditions  (fixed  grip  ox.  stress)  would  dissolve  some 
uncertainties. 
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ABSTRACT:  A  thermodynamic  description  of  fracture  based  on  percolation  theory  is 
studied  numerically.  The  fracture  behavior  under  fixcd-grip  and  constant  stress  conditions 
are  foimd  to  be  different.  Above  a  critical  fixed  strain  the  material  attains  an  equilibrium 
non-zero  crack  density  which  approaches  the  percolation  threshold  assymptorically  vdth 
increasing  strain.  This  equilibrium  may  or  may  not  be  readily  accessible  depending  upon 
the  applied  strain.  The  corresponding  fracture  stress  is  calculated  using  the  free  energy 
appropriate  to  constant  stress  conditions.  Effects  of  the  crack-crack  interaction,  anisotropic 
crack  density  and  pre-existing  pores  on  the  fracture  behavior  are  also  studied. 


I.  Ictrodcctxon 

Heccc5Jy,  zn  zUeazpi  hzs  bsen  szzie  [l^^]  to  rcfoinralate  tlie  GriSJli  critesxoa  for  cxzck 
groirth  \i]  fox  a  ^annlsx  rce^usa  c^g  i^ezs  Srosi  &e  ^emilailoa  tlicoiy  [5],  la  ilsis 
foranilaSioa  it  is  proposed  ilsat  t!:e  fosdiasnenfal  tbexiaodyaaaac'vanable  -splikli  serves  as 
&e  Older  parasaeter  is  tire  crack  coaceaSsalioa  p,  aad  aot  tke  leagtli  of  &e  sia^  or  &e 
doauaaat  cratk  as  in  tke  oiigbal  GriStli.  fireczj.  A  bask  cradr  has  a  ^icalsze  b,  which 
raaj  vary  betweeaiateratoauc  distances  to  raacroscopk  scales  (~  10~^ai)  depeadiag  apoa 
the  aratenaL  A  free  eaergy,  which  is  afaactioa  of  the  cradc  densify  onfy,  is  ihea  aimiaazed 
to  deteiaiiae  the  eqialibrhaa  crack  densify  ia  the  araterial.  Whea  a  material  of  voloase  V 
is  held  at  aa  average  straia  e,  (fixed  grip  conditions)  this  free  energy,  Fe{p)  caa  be  written 
ia  the  form 

J=i(p)  =  +u^+rfs-  Ts.  (1) 

Here  the  first  term  is  the  elastic  strain  energy  of  the  uncracked  material  whose  Young’s 
modulus  is  £<0,  U^m  is  the  mechanical  energy  released  under  fixed  grip  conditions  due 
to  the  reduction  of  the  stress  in  the  vicinity  of  the  crack,  Ug  is  the  excess  energy  in  the 
system  due  to  the  opening  of  new  crack  surfaces  and  S  is  the  configurational  entropy  of  the 
cracks.  If  on  the  other  hand,  the  external  stress  a  is  specified,  (constant  stress  conditions) 
the  appropriate  [6]  free  energy  that  should  be  minimized  is  the  Legendre  transformed  free 


energy 

FM=F,{p)-caV 

a^V  (2) 

=  -^  +  u^m  +  Us-ts, 

with  a  =  dF^ldi.  When  evaluated  on  a  lattice  representation  of  the  material,  with  a  lattice 
constant  6,  the  terms  appearing  in  the  free  energy  have  rather  simple  forms.  Denoting  the 
linear  size  of  the  material  by  L,  (in  units  of  6),  we  get  for  the  surface  energy 


Ug  =  2N^rp. 


2 


(3) 


Here  N  is  tlie  Bmc!>8r  of  feces  in  the  3  dimensional  lattice,  which  is  propoitipnal  to  and 
7  is  the  surface  energy  density  of  the  material.  In  a  2  dimensional  version  of  the  theory, 
N  would  be  the  number  of  bonds,  which  is  proportional  to  and  the  factor  I?  would  be 
replaced  hf  b.  For  the  case  of  externally  specified  stress,  the  mechanical  energy  released 


is  ^en  by  [7,8] 


while  for  a  system  with  specified  strain, 

(5) 

The  moduli  Ee  and  will  be  defined  later.  This  form  assumes  a  uniform  stress  a  (or 
strain  e)  throughout  the  material,  which  is  assumed  to  be  originating  from  uniaxial  tension. 
lij_  is  the  component  of  the  end-tc-end  distance  of  cluster  i  of  cracks  in  the  direction 
perpendicular  to  the  tension.  For  the  2  dimensional  case  the  exponent  3  has  to  be  replaced 
by  2.  The  magnitude  of  Uj^  is  a  monotonicly  increasing  function  of  p  up  to  the  percolation 
threshold,  in  the  vicinity  of  which  it  peaks  strongly.  This  pezdc  is  bounded  by  the  total 
voliune  (area  in  2  dimensions)  of  the  lattice  as  the  total  volume  within  which  stress  is 
relieved  cannot  be  greater  than  the  sample  volume. 

The  use  of  different  moduli  E^  and  E^  in  the  expressions  for  the  mechanical  energy 
released  is  essential  for  obtaining  the  same  effective  modulus  Egjf  under  both  conditions. 
As  one  can  define  Egfj  through  either 

=  (6) 


a'^V  o^V 


lEgfi  2Eo 

the  equivalence  of  both  definitions  prescribes  two  relations  between  the  three  moduli,  Ea^ 
Ee  and  Egj f.  We  are  thus  free  to  use  a  certain  physical  model  or  approximation  to  fix 
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one  of  those.  The  self  consistent  field  (SCF)  approximation  of  Budiansky  and  O’Connell 
[8]  results  in  the  choice 

Ea  =  Eo{l-v),  (8)  . 


where  v  =  follows  that  Ee  =  Eq  and  E^ff  =  E^. 

Another  choice,  which  can  be  shown  [3]  to  be  related  to  the  differential  scheme  (DS)  of 
Hashin  [9]  corresponds  to 


Ee  = 


Eq 

1  +  u’ 


(9) 


together  with  E^ff  =  Ee  and  Ea  =  Eq. 


A  parameterization  (which  is  not  unique)  that  includes  the  two  proceeding  results  as 
limiting  cases  is  to  write  down  the  three  moduli  in  terms  of  a  free  parameter  6  such  that 

Eq 


Ee  = 


1  +  dv  ’ 


E^  =  Eo[l-{l-e)v] 


(10) 


and 


Eeff  =  Eq 


1  -  (1  -  9)v 
1  +  Ov 


The  SCF  approximation,  Eq.  (8)  corresponds  to  5  =  0,  while  the  expression  resembling 
the  DS  approximation,  Eq.  (9)  is  equivalent  to  the  limiting  case  0  =  1. 


The  entropy  term  is  a  simple  combinatorical  contribution  resulting  from  the  number 
of  ways  of  placing  Np  cracks  on  N  faces.  It  is  given  by 


S  =  -Nks  (piogp  +  (1  -  p)  log(l  -  p)] .  (11) 


However,  since  k^T  <S  76^  except  for  atomic  scale  cracks,  the  entropy  term  can  be  ne¬ 
glected  for  all  practical  purposes. 


In  an  earlier  version  of  the  theory  (ij,  an  estimate  for  the  frzicture  stress  was  obtained 
using  some  simple  arguments.  Although  not  stated  explicitely  in  that  work,  the  form  of 
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the  free  energy  used  corresponds  to  that  of  DS  under  constant  stress  conditions.  This 
estimate  amounts  to  approximating  the  mechanical  energy  (Eq.  (4))  by  the  contribution 
from  the  largest  cluster  of  cracks  for  p  ~  pc-  When  one  plots  the  total  free  energy  as  a 
function  of  the  crack  density,  one  obtains  two  minima.  One  is  near  p  =  0  (for  kgT  <  76^) 
while  the  other  one  is  at  the  percolation  threshold,  pc-  At  low  values  of  the  stress,  the 
former  is  the  more  stable  one.  When  the  stress  is  larger  than  a  critical  ‘fracture’  stress, 
Of,  the  second  minima  becomes  more  stable.  Neglegting  the  contributions  to  from 
the  finite  clusters  of  cracks,  one  finds  that  <7/  =  \/6'yEoPc/6.  This  resembles  the  Griffith 
formula  (cTc  «  V'1^ol<^]  relating  the  critical  stress,  Cc  to  the  crack  length,  a;  however  the 
percolation  result  is  for  an  assembly  of  random  cracks  while  the  Griffith  theory  applies  to 
a  single  crack. 

There  are  several  points  in  this  theory  that  need  to  be  further  investigated.  First, 
approximating  Uj^  by  the  contribution  from  the  percolating  cluster  only  is  unfounded.  It 
may  be  the  case  that  the  stable  crack  density  is  obtained  at  a  value  of  p  different  from  the 
percolation  threshold,  in  which  case  one  needs  to  evaluate  the  sum  in  Eqs.  (4)  and  (5) 
rather  than  replacing  it  by  7,  its  upper  limit.  A  further  important  point  is  the  introduction 
of  crack-crack  interactions  into  the  theory.  This  point  has  been  answered  in  the  later 
version  of  the  theory  (2]  in  a  formal  manner.  However,  the  effects  of  the  interactions 
can  be  obtained  only  by  adetailed  calculation  involving  a  numerical  simulation.  Another 
important  generalization,  which  again  is  addressed  formally  in  Ref.  2  is  the  minimization 
of  the  free  energy  with  respect  to  a  position  and  orientation  dependen .  crack  density. 

In  this  work  we  address  these  points  using  numerical  simulations  of  two  dimensional 
version  of  the  model.  We  calculate  the  free  energies  defined  in  Eqs.  (l)  and  (2)  using 
definitions  of  the  moduli  Ef,  and  Ea  corresponding  to  both  SCF  and  DS  approximations. 
We  also  study  the  effects  of  interaction  between  cracks  by  adding  a  heuristic  term  to  the  free 
energy  to  represent  the  crack-crack  interaction.  We  also  partially  address  the  point  raised 
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in  the  preceeding  paragraph  about  the  anisotropic  crack  density  by  minimizing  the  free 
energy  with  respect  to  an  orientation  dependent  crack  probability,  which  assumes  different 
values  in  the  directions  parallel  and  perpendicular  to  that  of  the  uniaxial  tensile  stress.  In 
the  next  section  we  give  details  of  the  calculations,  including  the  form  of  the  interaction 
assumed.  Section  III  presents  the  results  for  the  cases  with  and  without  interaction.  Also 
presented  are  the  results  for  the  lattices  into  which  pre-existing  pores  have  been  introduced 
to  study  the  effect  of  pore  concentration  and  size  on  the  fracture  properties.  In  Section 
IV,  we  discuss  the  physical  relevance  of  these  results  and  compare  them  with  theory  where 
applicable.  The  results  are  summarized  briefly  in  Section  V. 

II.  Model  and  Method 

We  represent  the  material  in  two  dimensions  by  a  square  lattice  of  size  LxL.  The  cracks 
can  form  only  on  the  bonds  of  the  lattice.  The  probability  that  a  given  bond  represents 
a  crack  is  independent  of  the  status  of  other  bonds  and  is  given  by  p.  For  the  sake  of 
generality,  we  differentiate  between  bonds  that  are  parallel  to  the  direction  of  the  stress 
and  those  that  are  perpendicular  to  it  and  assign  them  crack  probabilities  pjj  and  p_L 
respectively. 

There  is  a  huge  number  of  realizations  of  the  lattice  consistent  with  a  given  (P||,Pjl)- 
Each  such  realization  corresponds  to  a  different  value  of  the  energy.  Since  the  only  order 
parameter  in  our  theory  is  the  pair  (p||,Px)>  each  realization  is  equally  allowed  in  the 
evaluation  of  the  average  free  energy.  (This  corresponds  to  an  infinite  temperature  since 
realizations  with  different  energies  have  equal  probabilities.)  To  calculate  this  average 
energy,  we  simulate  the  lattice  with  a  given  (pjj,pj.)  several  times  and  calculate  the  average 
free  energy  and  other  required  averages. 

The  simulations  are  performed  in  the  following  manner;  We  assign  to  each  bond  a 
random  number  between  0  and  1  and  place  a  crack  on  it  if  that  random  number  is  smaller 
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than  pj|  or  depending  upon  its  orientation.  We  then  identify  the  connected  cracks  that 
form  a  ‘cluster’  of  cracks.  Several  quantities  are  calculated  for  e^h  cluster;  these  are  the 
length  of  the  cluster  in  both  directions,  denoted  by  /||  and  l±,  coordinates  of  its  center 
of  mass,  its  radius  of  gyration  Rg  and  the  total  number  of  lattice  sites  it  covers.  These 
are  needed  to  calculate  the  different  terms  appearing  in  the  expression  for  the  free  energy 
corresponding  to  a  particular  crack  cluster.  J'ig.  1  shows  several  of  these  quantities  for  a 
typical  cluster. 

The  term  representing  the  excess  surface  energy  due  to  the  opening  of  cracks  is  cal¬ 
culated  trivially.  It  is  given  by  the  total  number  of  open  bonds  (cracks)  in  the  lattice, 
which  averages  to  +  Pj.)-  The  term  that  gives  the  amount  of  mechanical  energy 

released  needs  some  consideration,  since  some  of  the  cracks  fall  entirely  or  partly  inside 
the  region  from  which  stress  is  relieved  due  to  the  presence  of  another  crack.  Therefore 
simply  summing  over  of  each  cluster  of  cracks  will  overestimate  this  contribution.  To 
avoid  this  problem,  we  place  on  each  cluster  an  imaginary  square  of  side  centered  on 
the  center  of  mass  of  the  cluster.  A  counter  assigned  to  each  lattice  site  is  incremented  by 
1  if  the  site  is  inside  the  square,  by  0.5  if  it  is  exactly  on  an  edge  and  by  0.25  if  it  is  on  a 
corner  of  the  square.  A  site  may  be  included  in  more  than  one  such  square.  Any  counter 
that  reaches  the  value  1  is  not  changed  any  further.  We  then  sum  over  the  values  of  the 
counters  on  each  lattice  site  and  identify  this  sum  (after  multiplication  by  the  stress  or 
strain  dependent  factor)  as  the  strain  energy  released. 

The  elimination  of  the  overlap  of  the  stress-free  areas  around  each  crack  amounts  to 
a  certain  interaction  between  the  cracks;  it  is  a  short-range,  excluded  volume  type  of 
interaction.  However,  we  know  that  the  existence  of  a  crack  at  a  finite  distance  from 
another  one  changes  the  stress  intensity  at  the  tip  of  the  latter,  in  a  manner  dependent 
upon  their  length  and  relative  orientation  (10).  To  account  for  this  additional  long-range 
interaction,  we  assume  an  interaction  energy  between  any  pair  of  ‘crack  clusters’  in  the 
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lattice  with  a  magnitude  proportional  to 

=  (12) 

Here  is  a  vector  with  components  /y  and  of  the  craek  cluster  t,  /j  its  length,  is  the 

distance  between  the  centers  of  mass  of  the  two  clusters  and  a?  -  =  •.  The  form  of 

y  gi  93 

this  interaction  energy  is  rather  heuristic,  and  is  certainly  an  extreme  oversimplification 
of  the  reality.  Analytical  forms  for  crack-crack  interaction  are  known  only  for  some  very 
simple.crack  configurations  (10, 11).  For  the  composite  cracks  that  appear  in  our  percola¬ 
tion  formalism,  evaluation  of  an  interaction  energy  for  an  arbitrary  crack  configuration  is 
beyond  the  scope  of  any  realistic  calculation.  Another  important  effect  that  this  simple 
form  does  not  exhibit  is  the  screening  of  the  interaction  between  a  pair  of  cracks  by  a  third 
that  lies  on  the  line  connecting  the  centers  of  the  two.  We  have,  however,  chosen  to  use 
this  simple  form  with  the  following  considerations  in  mind:  The  magnitude  of  the  interac¬ 
tion  should  decrease  rapidly  with  the  distance  between  the  cracks,  while  the  length  scale 
characterizing  this  decay  should  depend  on  the  sizes  of  the  cracks.  The  interaction  is  also 
orientation-dependent  in  such  a  way  that  cracks  perpendicular  to  each  other  are  favored 
(negative  interaction  energy)  rather  than  cracks  with  parallel  orientation.  We  hope  that 
despite  its  obvious  drawbacks,  the  introduction  of  such  an  interaction  energy  will  capture 
the  basic  physics  of  interacting  crack  systems. 

In  a  given  configuration  of  the  lattice,  the  contributions  from  aii  pairs  are  summed 
over  and  the  result  is  averaged  over  several  configurations.  The  resulting  average,  which 
we  denote  by  J?/,  is  now  a  function  of  {p||,pj.)  only  and  it  modifies  the  ‘unperturbed’  free 
energy,  Fq  as 

F  =  Fo  +  A'^Ej  +  (13) 

that  is,  both  the  stress  dependent  and  stress  independent  terms  of  the  free  energy  are  mod¬ 
ified.  Here  F  and  Fq  may  be  for  fixed-grip  or  constant  stress  conditions..  The  assumption 
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of  the  same  form  of  the  interaction  for  the  stress-dependent  and  the  stress-free  interaction 
is  another  simplification  since  the  stress-free  interaction  may  be  shorter  ranged  than  the 
stress-dependent  one.  A  and  B  are  arbitrary  parameters  giving  the  relative  strength  of 
the  interaction  energy  compared  to  the  ‘unperturbed’  free  energy.  We  expect  that  this 
kind  of  long  range  and  size  dependent  interaction  will  not  contribute  significantly  to  the 
stress-independent  part  of  the  energy,  since  any  local  distortion  should  not  propagate  for 
more  than  several  atomic  distances  (I2j.  Thus  we  expect  A<.  B.  However,  we  keep  both 
contributions  and  calculate  the  interaction  energy  for  several  values  of  the  pair  A  and  B. 

in.  Results 

We  simulated  lattices  of  linear  size  L  =  20  and  50.  The  facts  that  the  the  number  of  crack 
pairs  that  need  to  be  considered  in  the  evaluation  of  the  interaction  term  grows  as  the 
square  of  the  immber  of  clusters  in  the  lattice  and  that  the  whole  simulation  has  to  be 
repeated  for  several  values  of  (p||,Px)  prevented  us  from  carrying  out  the  calculation  on 
larger  lattices.  It  seems,  however,  that  the  results  from  the  two  cases  do  not  indicate  an 
appreciable  finite  size  effect,  and  that  the  conclusions  can  be  safely  extrapolated  to  larger 
lattices. 

First,  let  us  present  the  results  for  the  case  in  which  p||  =  Pj.  =  p;  that  is,  bonds 
are  cracked  with  the  same  probability  regardless  of  their  orientation  with  respect  to  the 
direction  of  the  tensile  stress  acting  on  the  material.  Fig.  2  shows  the  efi'ective  Young’s 
modulus  (scaled  by  Eq)  as  a  function  of  p  for  the  two  approximations  corresponding  to 
0  =  0  and  0  =  1.  We  will  refer  to  the  two  cases  as  the  SCF  approximation  and  DS 
approximation  respectively,  although  these  names  refer  only  to  the  qualitative  resemblence 
between  our  two  limiting  cases  and  the  SCF  results  of  Budiansky  and  O’Connell  [8]  and 
the  DS  results  of  Hashin  (9).  The  SCF  effective  modulus  decays  linearly  for  p  <  1  and 
practically  vanishes  around  p  =  0.5.  A  quadratic  fit  in  the  range  0  <  p  <  0.2  yields 
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E.cff/Eo  =  1  —  0.98p  —  2.31p^,  which  if  extrapolated  vanishes  at  p  =  0.48,  very  close  to 
the  (bond)  percolation  threshold  of  the  lattice,  pc  =  0.5.  (In  the  simulations,  the  effective 
modulus  saturates  at  a  very  low  value  for  p  >  0.5  due  to  the  way  the  sum  in  Eqs.  (4)  and 
(5)  is  calculated.  As  the  center  of  the  spanning  cluster  of  cracks  does  not  always  coincide 
with  the  center  of  the  lattice,  the  imaginary  square  plciced  on  the  cluster  to  calculate  the 

4 

stress-free  area  does  not  cover  the  whole  lattice.  Thus  even  for  p  >  Pc  v/e  may  obtain 
t;  <  1,  leading  to  a  small  but  non-zero  effective  modulus.)  The  DS  modulus  starts  to  decay 
linearly  as  well,  but  saturates  at  Eo/2  for  p  >  0.5.  This  is  due  to  the  fact  that  v  will 
never  be  larger  than  1,  while  the  analogous  crack  density  parameter  of  Hashin  [9]  can  grow 
indefinitely. 

From  here  on  we  will  use  the  SCF  approximation  to  evaluate  the  free  energies,  unless 
otherwise  stated.  The  energies  are  measured  in  units  of  76  (cf.  Eq.  (3)).  Fig.  3  shows  the 
dependence  of  the  fixed-grip  free  energy  on  the  crack  concentration  for  several  values  of  the 
dimensionless  strain  variable  e  =  — without  crack-crack  interactions:  A  =  B  =^0.  The 
data  shown  are  for  L  =  50.  For  e  =  0,  the  free  energy  consists  only  of  the  excess  surface 
energy,  which  increases  linearly  with  the  crack  density.  As  e  is  increased,  a  local  minimum 
appears,  which  for  e  =  ei  a  1.0  (equivalent  to  e  =  ei  ci  y/A^lE^h)  attains  the  energy 
value  of  0  at  p  =  Pi  ~  0.40,  the  same  as  the  energy  value  in  the  intact  material  with  p  =  0. 
At  even  higher  values  of  the  strain,  this  local  minimum  becomes  the  absolute  minimum, 
which  gets  deeper  with  increasing  strain.  At  strains  for  which  e  >  62  2  (e  >  C2  —  V^ci)» 

the  initial  slope  of  the  free  energy  becomes  negative.  With  increasing  e  the  value  of  p  for 
which  Ff  is  minimum  moves  slightly  toward  higher  values  of  p,  approaching  the  critical 
value  of  Pc  =  1/2  for  much  higher  values  of  the  external  strain.  The  local  fluctuation  of  the 
curves  is  due  to  the  fact  that  only  20  realizations  of  the  lattice  were  simulated  for  each  p. 
A  similar  behavior  v/as  obtained  for  L  =  20,  where  400  lattices  were  used  for  the  averaging. 
The  results  for  ci  and  pi  were  slightly  different  from  the  ones  observed  for  L  —  50.  They 
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are  shown  in  Table  I  together  with  the  numerical  results  for  other  cases  analyzed. 

In  Fig.  4  we  show  the  average  constant-stress  free  energy  of  20  L  =  50  systems  as  a 
function  of  the  crack  concentration  at  several  values  of  the  applied  stress.  This  function 
exhibits  a  very  different  behavior  compared  to  that  of  the  fixed-grip  free  energy.  At  any 
finite  value  of  the  stress  the  free  energy  assumes  an  extremely  large  negative  value  as  p 
approaches  pc.  This  is  due  to  the  fact  that  within  the  SCF  approximation  the  effective 
modulus  vanishes  approciching  pc.  For  low  values  of  the  stress  the  energy  increases  at  low 
crack  concentrations  before  decrezising  rapidly  close  to  p^.  However  for  s  =  ^ 

{a  y/SEo'ilb)  this  positive  slope  vanishes  and  becomes  negative  for  higher  stresses. 

When  we  allow  the  crack  clusters  to  interact,  a  somewhat  different  picture  emerges. 
Fig.  5  shows  the  average  fixed  grip  free  energies  for  the  same  lattices  used  for  Fig.  3  with 
interaction  parameters  A  =  5  =■•  1.  For  e  =  0,  the  energy  again  rises  linearly  for  small 
values  of  p,  but  the  rise  becomes  much  steeper  when  the  interaction  between  the  cracks 
become  significant  at  larger  p.  It  seems  that  no  minimum  is  formed  now  at  a  non-zero 
value  of  p  before  e  c;  1.5;  rather  we  find  an  energy  that  increases  with  a  small  slope  with 
increasing  p,  before  it  shoots  up  at  around  p  a  0.3.  At  a  strain  that  corresponds  to  c  a:  1.5, 
a  global  minimum  is  observed  at  p  0.25.  It  seems  that  ei  and  e2  (corresponding  to  strains 
at  which  a  minimum  occurs  and  at  which  the  initial  slope  becomes  negative  respectively) 
coincide  in  this  case. 

We  calculated  the  free  energy  as  a  function  of  the  crack  concentration  for  several  values 
of  the  interaction  parameters  and  identified  ei  and  pi  for  each  case  as  described  above. 
An  increase  in  either  of  the  two  parameters  causes  ei  (and  consequently  the  corresponding 
strain)  to  increase  and  pi  to  decrease.  Table  I  summarizes  the  resulting  cj  and  pi  for  both 
L  =  20  and  L  =  50  as  a  function  of  the  interaction  parameters. 

The  restriction  of  py  =  p±  turns  out  an  unfavorable  one  for  our  system.  When  we  allow 
the  two  kinds  of  bonds  to  have  different  probabilities  to  represent  cracks,  the  energy  of  the 
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system,  is  lowered  by  decreasing  the  number  of  cracks  parallel  to  the  tensile  stress.  Li  fact, 
when  crack-crack  interaction  is  not  taken  into  account,  the  minimum  of  the  fixed-grip  free 
energy  at  any  given  stress  o  is  attained  on  the  p||  =  0  line.  Fig.  6  shows  this  behavior  for 
several  values  cf  c  on  the  line  pjj  =  0.  One  can  see  from  the  figure  that  the  local  minimum 
near  pj_i  ~  0.5  becomes  stable  at  e  =  ej  ci  0.58,  which  is  lower  than  the  one  obtained 
when  all  the  bonds  are  restricted  to  have 'the  same  probability  of  cracking.  This  means 
that  the  condition  represented  by  the  minimum  on  the  p|j  =  0  line  will  be  achieved  before 
any  other  local  minimum,  including  the  one  on  the  isotropic  line,  reaches  the  energy  value 
of  zero.  It  will  thus  be  the  relevant  one  whenever  the  condition  of  isotropicity  is  removed. 
As  in  the  isotropic  case,  at  even  higher  strains  the  Initial  slope  becomes  negative.  Here  we 
find  this  to  occur  for  c  >  63  0.8. 

One  would  expect  that  introduction  of  interactions  that  favor  cluster  pairs  that  are 
perpendicular  to  each  other  will  move  this  minimum  away  from  the  py  =  0  line.  We  find, 
however,  that  even  though  the  location  of  the  minimum  and  the  stress  level  at  which  it 
becomes  the  absolute  minimum  change,  it  still  remains  on  the  p||  =  0  line.  The  reason, 
apparently  is  that  when  we  introduce  a  small  number  of  cracks  parallel  to  the  axis  of 
the  tensile  stress  at  random,  some  of  them  go  to  regions  in  which  reduction  in  energy 
due  to  interactions  is  ineffective.  They  are  either  attached  to  already  existing  clusters 
of  ‘perpendicular’  cracks,  or  are  located  far  from  them,  so  that  the  magnitude  of  the 
interaction  is  small.  All  these  ‘parallel’  cracks,  however,  cost  additional  surface  energy  and 
the  overall  balance  is  unfavorable  to  the  system.  Only  for  values  of  p±  significantly  larger 
than  pj^i  do  we  see  a  favorable  effect  of  addition  of  ‘parallel’  bonds.  Using  A  =  B  =  1  for 
the  L  =  50  system,  we  find  that  an  absolute  minimum  is  obtained  at  p_L  2^  0.2  for  e  =  0.81, 
while  addition  of  ‘parallel’  bonds  starts  to  be  favorable  at  a  0.5. 

In  order  to  study  the  effects  of  porosity  on  the  fracture  behavior  within  our  model,  we 
repeated  our  simulations  with  lattices  on  which  pores  are  initially  placed.  A  pore  for  our 
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purpose  is  a  square  of  side  g  lattice  units  within  which  all  the  bonds  are  removed  (including 
those  on  its  perimeter).  Any  two  crack  clusters  that  touch  the  same  pore  are  assumed  to 
be  connected  for  the  purpose  of  calculating  the  strain  energy  released  and  the  interaction 
energy.  Any  crack  placed  at  random  inside  or  on  the  perimeter  of  the  pore  does  not  cause 
any  surface  energy.  In  accordance  with  Ref.  1  wc  studied  the  cases  with  g  =  1  (pore  size  ~ 
basic  crack  length)  and  g  =  5  (pore  size  >  basic  crack  length).  When  the  pore  size  is  equal 
to  the  crack  length,  we  find  that  for  pore  concentration  (f>  =  0.1,  the  absolute  minimum 
of  the  fixed-grip  free  energy  is  first  obtained  at  p  =  pi  ci  0.3  for  ci  oi  0.8  (with  isotropy 
restriction  and  in  the  absence  of  interactions).  pi  turns  out  to  be  smaller  than  the  one 
observed  for  the  (f>  =  0  system  since  the  crack  coalescence  is  assisted  by  the  pores.  62  is 
also  considerably  reduced  to  62  1.1.  The  introduction  of  interactions  again  changes  pi 

and  6i  in  a  manner  similar  to  the  pore  free  system.  Under  constant  stress  conditions,  we 
find  that  the  initial  slope  first  turns  negative  at  s  cs:  0.8.  Increasing  the  concentration  of 
the  pores  to  0.2,  we  obtain  ei  a  e2  ^  0.8,  pj  0.25  and  s  a  0.5.  When  the  pore  density 
of  0.1  is  achieved  with  pores  larger  than  the  basic  crack  size,  we  find  that  these  numbers 
approach  their  values  in  the  absence  of  any  pores.  For  pores  of  linear  size  of  5  lattice 
units,  we  find  that  ei  cs  1.0.  pi  ci  0.4,  62  1.5  and  5  ~  1.2,  approaching  the  corresponding 

values  for  the  pore-free  lattice.  Thus  our  data  confirm  the  prediction  made  in  Ref.  1. 


13 


47 


IV.  Discussion 

One  interesting  result  of  the  simulations  is  the  appearance  of  the  minimum  of  the  fixed- 
grip  free  energy  at  a  crack  concentration  significantly  below  the  critical  concentration  for 
percolation  when  a  self  consistent  field  approximation  is  used.  This  is  true  both  with  and 
without  the  interaction  terms  and  remains  such  when  the  isotropic  restriction  is  removed, 

4 

since  with  pj|  =  0,  the  critical  value  for  is  1.  The  implication  is  that  at  strains  c  that 
correspond  to  the  dimensionless  strain  variable  ej  <  c  <  C2>  there  exists  a  stable  crack 
concentration  for  the  material.  However,  if  one  starts  with  initially  crack-free  material, 
this  minimum  is  inaccessible  since  the  free  energy  is  an  increasing  function  of  the  crack 
concentration  at  that  strain  range.  Only  when  e  >  62  can  this  minimum  be  reached 
spontaneously.  Fracture,  that  corresponds  to  an  equilibrium  crack  density  of  Pc  can  then 
occur  only  at  an  infinite  strain  since  the  minimum  approaches  Pc  only  asymptotically  with 
increasing  e. 

This  result  is  in  fact  not  surprising  since  in  our  model  fracture  is  due  to  the  coalescence 
of  random  crack  clusters,  which  can  occur  only  at  p  =  pc.  Because  within  the  self-consistent 
field  approximation  the  effective  Young’s  modulus  vanishes  at  that  concentration,  any  finite 
fracture  stress  would  correspond  to  infinite  strain.  In  order  to  know  the  fracture  stress, 
one  needs  to  analyze  the  system  under  constant-stress  conditions.  From  Fig.  4  we  find  that 
indeed  such  a  finite  fracture  stress  exists.  In  our  model  system,  at  stresses  that  correspond 
to  the  dimensionless  stress  variable  s  >  2.0  the  material  would  spontaneously  break.  Thus 
the  fracture  stress  in  our  model  is  given  by  ay  oi  \/8')EoJb.  At  this  stress,  the  intact 
material  has  an  elongation  that  corresponds  to  63  =  cry/j&o  since  the  effective  Young’s 
modulus  for  the  intact  material  is  Eq.  If  kept  at  this  constant  stress,  the  equilibrium  state 
of  the  material,  which  is  at  p  =  pc  would  be  reached  at  a  strain  that  increases  without 
bound.  Thus  in  a  sense  our  model  describes  a  perfectly  plastic  fracture. 

Within  the  SCF  approximation,  the  Young’s  modulus  vanishes  linearly  at  the  percola- 
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tioa  thiesbo!d.  The  apgumfeaisaa  ihzi  lesez^les  SEe  DS  asgjcem  ^es  a  fe!te  raize  kk 
tEe  elastic  cozstant  at  tliegeoiisetricpegcoIaticntEresEoId.  Ib  f2dt.3  MasEEi’s  feeataaezt  for 
a  cracked  materia!  [S^  does  cot  take  mto  cczrideratircr  a  Scife  peecola^mi  tEresboId  ^ 
cracks  of  c^Hg^Ie  voEuoe  acd  piedkts  ir^oG  feat  ranisE  at  an  iaaaite  cozcentrar 

tion  of  such  cracks.  Treatmezts  feat  go  bejozd  fee  eSective  ssedioaa  appsorsmaiioa  (wMch 
both  approaches  utilize)  g^e  a  soft  decay  of  fee  modalzs  at  fee  percolation  feresEold  pc 
of  fee  materia!  as  E^fj  cc  {p—pcY ,  wife  fee  best  estimate  [13,14]  for  T  beaug  Tc^4.  bofe 
in  2  and  3  dimensions.  (Li  terms  of  fee  crack  cozceniration,  fee  material  fragmentation 
threshold  pc  [15,16]  corresponds  to  what  has  been  termed  [l7j  p[  or  p^jj  of  the  csarls, 
in  contrast  to  their  percolation  threshold  termed  for  fee  present  purposes  p^j  (<  ?cjj,  la 
three  dimensions) 

It  is  interesting  to  note  that  a  recent  study  by  Blamberg-Selinger  ct  aJ  [l8]_feat 
utilized  a  thermodynamic  approach  to  fracture  also  resulted  in  an  equilibrium  state  with  a 
finite  density  of  defects.  That  study  involved  the  relaxation  of  an  elastic  network  in  which 
equilibrium  defects  are  thermally  formed  or  destroyed  until  an  equilibrium  configuration 
is  reached.  At  a  certain  range  of  the  applied  stress  or  strain,  the  equilibrium  configuration 
is  one  in  which  a  finite  fraction  of  the  lattice  bonds  are  broken.  This  equilibrium  may 
be  stable  or  metastable.  In  that  study,  however,  the  fracture  is  essentially  a  brittle  one, 
unlike  our  model  in  which  it  results  from  the  coalescence  of  random  crack  clusters 

The  introduction  of  a  long  range  interaction  between  cracks  changes  the  fracture  be¬ 
havior  qualitatively.  It  seems  that  under  fixed-grip  conditions,  the  distinction  between 
and  C2  vanishes.  This  means  that  the  there  is  no  local  minimum  of  the  free  energy  at  non¬ 
zero  values  of  the  crack  concentration.  We  also  find  that  irrespective  of  the  parameters  of 
the  interaction  strength,  the  introduction  of  crack-crack  interactions  increases  the  value  of 
Cl.  Even  though  the  interaction  between  any  two  cracks  can  have  a  positive  or  negative 
effect  on  the  energy  of  the  system,  the  averaging  process  seems  always  to  result  in  the 
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same  trezd  for  e«  or  for  the  &acteie  stress. 

An  jaJeresfeg  resaS  of  fire  modd  is  the  instzbilsfy  with  respect  to  the  zssymetiy  in 
the  piohabahires  of  cradc  formation  in  the  directions  parallel  and  perpendicular  to  the 
direction  of  the  external  stress.  The  fact  that  in  all  cases  studied  fracture  occurs  when 
the  densify  of  cracks  in  the  direction  parallel  to  the  stress  is  zero  means  that  these  cracks 
cost  energy  to  form  and  do  not  relieve  considerable  stress  energy,  even  when  interactions 
are  taken  into  acconnt.  This  behavior  is  the  closest  that  one  can  approach  brittle  macture 
vnthin  our  model,  which,  due  to  its  probabilistic  nature,  does  not  acconnt  for  the  growth 
of  a  single  crack.  In  realify  we  would  expect  that  the  microscopic  structure  of  the  material 
vrill  allow  for  some  relation  or  restriction  between  the  densities  in  different  directions  in 
space  (vTcbably  as  a  continous  function  of  the  direction).  Li  our  discrete  model,  this  would 
correspond  to  looking  for  the  minimum  of  the  free  energy  along  some  line  intermediate 
between  the  isotropic  line  and  the  pg  =  0  line. 

The  dependence  of  the  fracture  stress  on  the  pore  density  and  pore  size  is  a  feature 
predicted  by  the  theoretical  analysis  of  the  model  (ij.  The  behavior  found  in  this  work 
is  in  qualitative  agreement  with  that  prediction,  although  the  present  study  relies  on  the 
SCF  approximation  while  Ref.  1  used  the  DS  approximation.  According  to  the  theory, 
when  the  pore  size  is  of  the  same  order  of  magnitude  as  the  basic  crack  size,  the  fracture 
stress  should  decrease  as  oy  oc  -y/(l  —  ^)Pc(^).  ^  is  the  density  of  pores  and  pc(^)  is  the 
percolation  threshold  of  the  cracks  in  a  lattice  with  a  density  (j>  of  pores-  Li  our  simulations 
we  found  the  fracture  stress  to  decrease  froin  what  corresponds  to  s  =  2  for.  a  pore  free 
system  to  s  =  0.8  and  0.5  for  systems  with  ^  =  0.1  and  0.2  respectively.  We  also  verify 
numerically  the  prediction  that  as  the  pore  size  increases,  the  fracture  stress  increases 
towcirds  its  value  in  pore-free  material.  For  ^  =  0.1  increasing  the  pore  size  to  5  lattice 
units  was  found  to  increase  the  frticture  stress  to  1.2.  This  behavior  is  also  experimentally 
observed  in  several  ceramic  materials  (l9j. 
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V.  Summary 

A  thermodynanuc  description  of  fracture  due  to  the  coalescence  of  cracks  and  pores  seems 
to  show  several  interesting  features  when  investigated  using  percolation  techniques.  We 
find  that  one  obtains  different  results  if  the  fracture  is  attained  under  fixed-grip  or  constant- 
stress  conditions.  When  the  material  is  held  at  a  fixed  strain  a  stable  configuration  with 
a  crack  density  pi  oi  0.4  is  formed  for  strains  above  cj,  which  in  a  pore  free  material  with 
no  crack-crack  interaction  is  given  by  =  ■\/A'^IEqo.  This  state  becomes  spontaneously 
accessible  at  a  strain  £3  =  V2£l-  The  corresponding  crack  density  increases  from  pi  to 
Pc  asymptotically  with  increasing  strain.  Fracture  occurs  only  at  infinite  strain.  The 
fracture  stress  is  finite,  however,  since  the  effective  Young’s  modulus  Veinishes  linearly 
approaching  the  percolation  threshold.  This  stress,  oj,  can  be  found  using  the  free  energy 
appropriate  for  the  constant  stress  conditions.  We  find  that  it  is  given  by  try  =  y/8'^Eofb  = 
£^0^2*  Crack-crack  interaction  as  introduced  using  a  heuristic  form  causes  the  distiction 
between  ej  and  £2  to  disappear,  thereby  rendering  the  equilibrium  crack  concentration 
spontaneously  accessible. 

When  crack  densities  parallel  to  the  direction  of  the  tensile  strain  and  perpendicular 
to  it  are  allowed  to  be’  different,  the  system  achieves  its  minimum  free  energy  by  having 
all  the  cracks  in  the  direction  perpendicular  to  the  direction  of  the  strain.  Furthermore, 
this  minimum  is  attained  at  a  smaller  applied  strain  than  the  one  needed  to  achieve  the 
metastable  state  with  isotropic  crack  densities.  Introduction  of  the  crack-crack  interaction 
does  not  alter  this  situation  cither. 

For  systems  that  contain  pores,  the  fracture  energy  is  considerably  reduced.  The 
fracture  stress  is  also  found  to  be  an  increasing  function  of  the  ratio  of  the  pore  size  to  the 
basic  crack  length,  in  agreement  with  the  situation  observed  in  certain  ceramic  materials. 
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Table  I 

ei  and  pi  as  a  function  of  the  crack-crack  interaction  parameters  A  and  B  for  lattices  of 
size  L  =  20  and  50,  The  dimensionless  strain  variable  ei  =  gives  the  strain  ei  at 
which  the  local  minimum  in  the  fixed-grip  free  energy  at  a  nonzero  p  becomes  the  absolute 
minimum.  pi  is  the  corresponding  crack  concentration, 

L  =  20  L  =  50 


A 

B 

ei 

Pi 

ei 

Pi 

0.0 

0.0 

1.13 

0.41 

1.00 

0.40 

0.0 

0.5 

1.32 

0.35 

1.20 

0.34 

0.0 

2.0 

1.68 

0.22-0.24 

1.57 

0.22 

0.0 

5.0 

1.96 

0.01-0.10 

1.82 

0.11 

0.5 

0.0 

1.29 

0.37 

1.17 

0.34 

2.0 

0.0 

1.54 

0.26-0.28 

1.44 

0.28 

5.0 

0.0 

1.80 

0.16-0.19 

1.70 

0.14 

0.5 

0.5 

1.43 

0.28-0.31 

1.31 

0.28 
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FIGURE  CAPTIONS 


Figure  1.  Definition  of  several  of  the  quantities  calculated  for  a  typical  crack  cluster 
with  11  bonds  (cracks),  x  shows  the  position  of  the  center  of  mass  of  the  cluster.  The 
radius  of  gyration  of  this  cluster  is  1.87  lattice  units,  a  is  the  external  tensile  stress. 

Figure  2.  The  effective  Young’s  modulus  at  crack  concentration  p  within  (a)  self 
consistent  field  approximation  and  (b)  differential  scheme. 

Figure  3.  Average  fixed-grip  free  energy  of  20  lattices  of  linear  size  L  =  50  as  a  function 
of  the  crack  density  for  several  values  of  the  dimensionless  strain  variable  e  = 
indicated  next  to  each  curve.  The  interaction  parameters  A  =  B  =  0. 

Figure  4.  Average  constant-stress  free  energy  of  20  lattices  of  linear  size  L  =  50  as  a 
function  of  the  crack  density  for  several  values  of  the  dimensionless  stress  variable  s  = 
indicated  next  to  each  curve.  The  interaction  parameters  A  =  B  =  0. 

Figure  5.  Average  fixed-grip  free  energy  of  a  lattice  of  linear  size  L  =  50  as  a  function  of 
the  crcick  density  for  several  values  of  the  dimensionless  strain  variable  e  =  indicated 
next  to  each  curve.  The  interaction  parameters  A  =  B  =  1. 

Figure  6.  Average  fixed-grip  free  energy  versus  density  of  cracks  in  the  direction 
perpendicular  to  the  axis  of  the  tensile  stress.  The  density  of  cracks  in  the  other  direction 

p  I 

is  zero.  The  value  of  the  dimensionless  strain  variable  e  =  is  shown  next  to  each 
curve. 
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ABSTRACT: 

We  propcse  several  models  for  crack  propagation  in  ceramic  materials  with  microcrack 
formation  and  coalescence.  The  models  differ  in  the  relative  time  scales  involved  in  the 
internal  relaxation  processes  and  the  external  loading  rate.  We  find  that  when  the  external 
loading  rate  is  very  high,  a  rising  R-curve  is  obtained.  For  slow  external  loading,  the  results 
depend  upon  the  relative  rate  of  the  different  internal  relaxation  processes  and  upon  the- 
distribution  of  the  critical  stresses  needed  to  form  the  micrccracks.  In  this  regime,  the 
existence  of  a  rising  R-curve  behavior  is  found  to  be  primarily  due  to  microcrack  shielding. 


L  INTRODUCTION 
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The  strength  of  brittle  materials  is  limited  by  the  existence  of  inherent  flaws  within  the 
materials.  A  quantitative  statement  of  this  observation  was  first  given  by  Griffith^,  who 
coiisidered  an  infinite  plate  of  uniform  thickness  in  which  a  thin  crack  of  length  2a  is 
inserted.  Using  a  stress  analysis  due  to  Inglis^,  GriflSth  showed  that  when  the  material  is 
under  a  tensile  stress  a  applied  at  infinity,  an  increase  of  Ac  in  the  crack  length  will  reduce 
the  strain  energy  of  the  material  by  2w6c^aAc/5,  where  E  is  Young’s  modulus  and  b  the 
thickness.  However,  the  increase  in  the  crack  length  will  cost  a  surface  energy  of  26Aa7. 
Here  7  is  the  surface  energy  density  of  the  solid  and  the  factor  2  comes  from  the  formation 
of  two  new  surfaces.  The  crack  length  will  spontaneously  increase  when  the  reduction  in 
the  strain  energy  becomes  larger  than  the  increase  in  the  surface  energy.  This  gives  the 
famous  Griffith  criterion  for  the  fracture  stress,  a j: 

Of  =  \/2E'tla-n  .  (1) 

One  can  interpret  R  =  27  as  the  material’s  resistance  to  further  fracture.  Experimental 
test  of  this  prediction,  performed  by  Grilfith  himself  on  glass,  verified  the  dependence  on 
the  flaw  size,  but  the  critical  stresses  turned  out  to  be  several  times  higher  than  the  pre¬ 
dicted  value.  The  deviation  was  found  to  be  larger  for  metals^.  The  discrepancy  implies 
that  there  is  another  sink  for  energy  in  addition  to  the  new  crack  surfaces.  Irwin^  sug¬ 
gested  that  this  sink,  at  least  for  metals,  is  a  plastically  strained  zone  around  the  crack 
tip  within  which  energy  is  dissipated.  For  such  cases,  the  resistance  to  frcicture,  R  is 
composed  of  the  sum  of  the  surface  energy  and  the  plastic  energy.  In  cerzimic  materials 
and  some  multiphase  materials,  like  concrete,  in  which  plasticity  is  not  significant,  there 
are  additional  mechanisms  that  dissipate  energy  during  fracture.  A  prominent  toughening 
mechanism,  observed  in  zirconia  and  alumina  ceramics,  is  the  so-called  ‘transformation 
toughening’^.  Here  energy  is  dissipated  in  the  transformation  of  the  tetragonal  Zr02  par¬ 
ticles  dispersed  in  the  matrix  into  a  monoclinic  structure.  Another  mechanism  suggested 
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for  non-transforming  ceramics  is  due  to  a  restraining  force  across  the  newly  formed  crack 
interface  caused  by  grain  bridging®*®.  Quoting  some  earlier  works^*®*®,  the  authors  sug¬ 
gest  that  the  same  mechanism  may  be  relevant  to  other  materials  like  steel,  concrete  and 
rocks.  A  third  important  mechanism,  on  which  we  will  base  our  study,  is  the  microcrack 
toughening.  Here  a  cloud  of  microcracks  that  form  in  front  and  at  the  sides  of  the  main 
crack  absorb  some  of  the  energy  that  would  otherwise  be  used  in  increasing  the  length  of 
the  macrocraick  in  a  brittle  fracture.  They  may  also  act  as  shields  that  screen  the  macro¬ 
crack  tip,  thereby  reducing  the  stress  intensity  there.  Substantial  experimental  evidence 
to  the  relevance  of  microcracking  mechanism  was  provided  by  the  works  of  Knehens  and 
Steinbrech^®  and  Ruble  et  ai 

These  mechanisms  have  also  been  used  to  explain  an  interesting  feature  observed  in 
ceramic  materials,  the  so-called  R-curve  behavior.  It  is  found  in  these  materials,  that 
when  one  evaluates  the  resistance  to  fracture  R  defined  above  as  a  function  of  the  crack 
length  increment  by  some  standard  procedures^^*^^,  one  obtains  a  resistance  that  increases 
with  increasing  crack  length.  Such  a  behavior  has  immense  engineering  significance,  as  the 
implication  of  this  increase  is  that  unstable  propagation  of  the  crack  is  delayed  until  the 
crack  reaches  larger  sizes,  rendering  the  detection  of  the  crack  before  a  catastrophic  failure 
possible.  Theoretical  models  based  upon  the  previously  mentioned  mechanisms®*  have 
concluded  that  a  rising  R-curve  is  possible  through  these  mechanisms.  A  model  treating 
the  formation  of  individual  microcracks  was  proposed  by  Hoagland  and  Embury^®  who 
simulated  the  microcracking  process  around  a  long  crack  by  solving  for  the  stresses  in 
the  material  in  a  self-consistent  manner.  They  found  the  ninnber  of  microcracks  formed 
for  several  possible  distributions  of  microcrack  formation  stresses  and  orientations  as  a 
function  of  external  stress  intensity.  R-curves  for  these  distributions  were  also  calculated 
under  the  assumption  that  the  material  still  exhibits  linear  elastic  behavior.  With  this 
assumption,  the  fracture  resistance  is  proportional  to  where  Kj  is  the  stress  intensity 
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factor®.  The  R-curves  turn  out  to  be  continously  increasing.  However,  due  to  the  intensive 
calculations  involved  in  the  determination  of  the  stress  distribution,  only  crack  increments 
of  about  12  microcrack  length  units  could  be  realized.  The  model  was  later  extended^^  to 
cases  where  both  microcracking  and  stable  crack  extension  occur  simultaneously. 

In  this  work  we  suggest  three  models  similar  to  the  Hoagland-Embury  model.  We 
believe  that  each  of  these  models  is  applicable  to  cases  with  different  relative  rates  of 
the  internal  relaxation  processes  of  the  material  and  of  loading.  To  our  knowledge,  ear¬ 
lier  studies  of  fracture  or  other  breakdown  phenomena  do  not  take  these  different  rates 
into  account^®.  We  calculate  the  stresses  in  a  very  simplified  manner,  allowing  for  the 
cancellation  of  the  crack-microcrack  interaction.  This  makes  it  possible  to  separate  the 
contribution  to  the  R-curve  behavior  from  different  sources:  energy  dissipation  due  to 
microcrack  formation  and  microcrack  shielding.  The  R-curve  is  also  evaluated  using  two 
different  methods  as  these  may  lead  to  slightly  different  R-curve  behavior^®’^®.  In  the  next 
section  we  define  the  different  models,  stating  the  rate  regimes  for  which  we  believe  they 
are  relevant.  In  section  III,  the  results  for  the  different  models  with  several  random  distri¬ 
butions  for  the  microcrack  formation  stresses  are  presented.  We  discuss  the  significance  of 
the  results  in  section  IV,  and  identify  the  conditions  under  which  a  rising  R-curve  should 
be  expected.  We  then  give  a  brief  summary  of  the  study  in  the  final  section. 
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n.  MODEL 


We  begin  with  a  two  dimensional  intact  material  lying  in  the  xy-plane  in  which  a  semi¬ 
infinite  crack  on  the  negative  x-axis  is  initially  inserted.  The  material  contains  regions 
in  which  a  microcrack  may  potentially  be  formed.  For  simplicity,  we  assume  that  these 
potential  microcracks  are  located  in  a  regular  manner  as  the  bonds  of  a  square  lattice. 
The  length  of  a  bond  defines  the  length  scale,  which  is  equivalent  to  the  typical  micro- 
crack  size.  Our  simplification  amounts  to  neglecting  the  variations  in  the  length  of  the 
microcracks  that  may  be  formed  and  to  assiuning  only  two  possible  directions  for  the 
microcracks.  Furthermore,  since  two  neighboring  microcracks  are  in  contact,  the  model 
allows  for  coalescence  of  individual  microcracks  among  themselves  or  with  the  main  crack. 


We  model  the  evolution  of  the  cracking  process  in  the  material  in  the  following  manner: 
The  material  is  placed  under  a  uniaxial  tension  in  the  y-direction  applied  at  infinity.  The 
stress  distribution  in  the  material  is  then  mainly  governed  by  the  singularity  at  the  tip  of 
the  long  crack.  For  this  so-called  Mode  I,  the  stresses  at  any  point  in  the  vicinity  of  the 
crack  is  given  by^^ 


A  ^  5  .  30 

1  +  sin  -sm  — 

V  2  2 


) 


/  Of  .  0  .  zo\ 

Here  (r,0)  are  the  polar  coordinates  of  a  point  in  a  coordinate  system  centered  on  the 
crack  tip.  The  stress  intensity  Kj  is  related  to  the  external  stress  a°°  and  to  the  creick 
length  a  as  Kj  =  a°°y/a.  If  microcracking  has  already  occurred,  the  stress  intensity  may 
be  modified  due  to  the  screening  by  the  microcracks,  as  will  be  explained  below.  Each 
bond  in  the  lattice  has  a  critical  stress  such  that  when  the  stress  on  the  bond  exceeds 
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that  critical  stress,  a  microcrack  is  irreversibly  formed.  The  critical  stresses  are  chosen 
randomly  from  a  given  distribution.  Initially  we  start  with  a  small  external  stress  and 
increase  it  until  microcracking  starts  to  occur.  The  exeu:t  rules  by  which  the  microcracks 
form  will  be  given  below.  During  the  course  of  the  simulation,  some  of  the  microcracks 
join  the  main  crack,  causing  it  to  extend  to  positions  with  positive  x.  We  monitor  the 
instantaneous  position  of  the  crack  tip,  given  by  the  point  with  the  largest  value  of  the 
x-coordinate.  That  point  is  used  as  the  origin  of  the  coordinate  system  in  the  calculation 
of  the  new  stress  distribution.  However,  even  when  that  tip  position  is  not  on  the  x-axis, 
the  overall  direction  of  the  crack  is  still  taken  to  be  horizontal,  since  the  initial  crack  was 
assumed  to  be  semi-infinite  in  that  direction.  This  assumption  results  in  the  crack  growing 
roughly  in  the  horizontal  direction,  with  slight  fluctuations  in  the  vertical. 

At  a  given  stage  of  the  simulation,  the  material  contains  several  isolated  microcracks 
and  clusters  of  microcracks  that  are  not  linked  to  the  main  crack.  Their  presence  alters  the 
stress  distribution  in  the  material.  In  particular,  the  stress  intensity  at  the  tip  of  the  main 
crcick  is  changed.  Kachanov^^  has  found  that  the  interaction  between  a  macrocrack  and 
one  or  more  microcracks  is  extremely  dependent,  both  qualitatively  and  quantitatively, 
on  their  relative  configuration.  The  stress  intensity  at  the  crack  tip  may  be  enhanced  or 
reduced  depending  upon  the  configuration.  However,  it  is  found  that  reduction  is  more 
dominant,  due  to  the  shielding  of  the  external  stress  by  the  microcracks.  He  found  that 
with  a  pair  of  microcracks  that  htve  the  same  orientation  as  the  main  crack  and  are 
symetrically  located  above  and  below  the  tip  of  the  horizontal  crack  (with  the  direction 
of  the  uniaxial  stress  being  vertical)  the  stress  intensity  at  the  crack  tip  is  shielded.  The 
stress  intensity  vanishes  rapidly  with  decreasing  distance  between  the  main  crack  and  the 
microcracks  (Fig.  3b  in  Ref.  22).  We  take  crack-microcrack  interaction  into  ciccount  only 
for  this  type  of  configuration.  In  our  lattice  representation,  this  includes  all  the  horizontal 
microcracks  whose  either  tip  has  the  same  x-ccordinate  as  the  macrocracL  Since  no  closed 
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form  for  the  dependence  of  this  interaction  on  the  crack-microcrack  separation  is  known, 
we  assume  that  the  reduction  in  the  stress  intensity  fcictor  at  the  tip  of  the  macrocrack  is 
given  by 

^  [viic-ytipf' 

y^c  and  yup  are  the  y-coordinates  of  the  microcrack  labeled  as  fic  and  of  the  crack  tip 
respectively.  The  sum  is  only  on  those  microcrzicks  that  are  in  the  configuration  described 
above,  a  is  a  parameter  that  can  be  varied  between  0  and  1  to  control  the  relative 
importance  of  the  shielding.  The  factor  S/tt^  assures  that  for  a  =  1  the  stress  intensity 
at  the  tip  is  completely  shielded  when  an  infinite  stack  of  microcracks  are  formed  above 
and  below  the  tip.  The  additive  property  of  Eq.  (3)  is  certainly  an  oversimplification, 
but  we  hope  that  this  simple  form  will  capture  the  essential  physics  involved  in  the  crack- 
microcrack  interaction  regarding  the  process  of  microcrack  toughening. 

Several  models  differening  in  the  exact  breaking  rules  of  the  potential  microcracks  are 
studied.  Physically,  these  models  differ  in  the  relative  magnitudes  of  the  time  scales  in¬ 
volved  in  the  three  processes  taking  place  simultaneously  in  the  system:  time  of  formation 
of  microcracks  in  the  overstressed  regions  (r^c)>  relaxation  time  of  the  local  stress  distribu¬ 
tion  when  the  external  stress  is  changed  or  when  the  crack  tip  moves  (r/)  and  the  inverse 
rate  of  increase  of  the  external  stress  (r^). 

MODEL  A:  For  a  given  external  stress,  ct°°,  calculate  the  stresses  on  all  the  bonds  on  a 
lattice  of  linear  size  L.  If  <7yy  on  a  horizontal  bond,  or  axx  on  a  vertical  one  is  above  the 
critical  stress  assigned  to  that  bond,  the  bond  is  broken  (a  microcrack  is  formed  there). 
The  connectivities  of  the  microcracks  are  established  and  the  new  crack  tip  is  identified. 
The  tip  is  the  point  on  the  main  crack  with  the  highesc  value  of  x-coordinate.  The  stress 
level  is  increased  by  a  prespecified  amount  and  the  stress  distribution  in  the  material  is 
calculated  again  with  the  new  tip  position.  The  process  is  continued  until  the  crack  tip 
reaches  a  boundary  of  the  lattice.  We  think  that  this  model  is  appropriate  for  the  case  of 
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very  fast  loading  rates  in  which  re  <  r^c  <  r/.  Although  for  most  brittle  elastic  materi¬ 
als  at  room  temperature  these  conditions  are  not  satisfied,  viscoelastic  materials  such  as 
polymeric  materials  or  metals  and  ceramics  at  high  temperatures  are  possible  candidates 
for  this  model^^. 

MODEL  B:  This  is  similar  to  Model  A  except  that  if  the  tip  position  changes  beacuse  of 
microcrack  coalescence,  stress  distribution  is  calculated  again  without  raising  the  external 
stress  level.  This  is  repeated  until  the  position  of  the  crack  tip  does  not  change  any  further. 
Only  then  is  the  external  stress  level  increased.  This  model  is  therefore  appropriate  when 
the  local  stresses  relax  more  rapidly  than  the  rate  of  increase  of  the  external  stress  but 
more  slowly  than  the  formation  of  a  microcrack  (r^ic  <ti  <  Te). 

MODEL  C;  After  identifying  the  bonds  on  which  the  stress  is  above  the  critical  stress  (as 
in  Models  A  and  B)  only  the  bond  with  the  largest  excess  stress  is  broken.  If  this  does  not 
change  the  crack  tip  position,  the  bond  with  the  next  largest  excess  stress  is  broken.  If  a 
microcrack  formed  changes  the  position  of  the  tip,  the  new  stress  distribution  is  calculated 
and  the  bond  with  the  largest  excess  stress  is  broken.  This  model  is  relevant  when  the 
relaxation  time  of  the  local  stresses  is  shorter  than  the  time  of  formation  of  a  microcrack, 
which  in  turn  is  shorter  than  inverse  rate  of  increase  of  the  external  stress,  r/  <  <  ^e- 


8 


71 


m,  RESULTS 

During  the  course  of  the  simulations,  we  monitor  the  number  of  microcracks  formed  {N) 
cuid  the  increase  in  the  length  of  the  macrocrack  (Ao)  as  a  function  of  the  external  stress 
(or  equivalently  of  the  nominal  stress  intensity  factor,  Kj).  Using  these  data,  we  estimate 
the  resistance  to  fracture  using  two  different  commonly  accepted^^’^^  measures  of  that 
property.  The  first  measure  is  the  rate  of  change  of  the  excess  surface  energy  of  the 
system  with  increasing  crack  length.  In  our  system,  the  excess  surface  energy  is  due  to 
the  microcracks  formed,  so  that  the  resistance,  R  is  given  by  R  =  2'ydNf  da.  Here  7  is  the 
surface  energy  density.  For  ideally  brittle  materials,  the  new  surfaces  formed  are  those  at 
the  tip  of  the  crack  and  this  definition  leads  to  a  constant  value  o{  R  —  27,  independent 
of  crack- length.  The  other  measure  of  R  is  based  upon  the  assumption  of  linear  elastic 
behavior  in  the  material.  In  this  case,  R  oc  ifj,  with  the  proportionality  constant  being 
a  geometry  dependent  combination  of  the  Young’s  modulus  and  Poisson’s  ratio^^.  These 
two  measures  are  equivalent  for  brittle  fracture;  however,  some  toughening  mechanisms 
are  found  to  lead  to  qualitatively  different  behavior  for  these  two  definitions^®.  We  expect 
similar  discrepancies  with  the  microcracking  mechanism  and  thus  plot  both  measures  of 
resistance. 

Fig.  1  shows  examples  of  the  crack  structures  generated  by  the  three  models  when  the 
macrocrack  traversed  a  square  lattice  of  linear  size  60  lattice  units.  The  initial  tip  of  the 
macrocrack  is  on  the  leftmost  bond  of  the  lattice.  In  all  the  three  examples  shown,  the 
effects  of  the  microcrack  shielding  is  ignored  (a  =  0).  The  critical  stresses  are  taken  from 
a  Weibull  distribution^^, 

<  cr)  =  1  -  exp  {-(a/oo)'”)  ,  (4) 

where  the  left  hand  side  denotes  the  probability  that  the  critical  stress  of  a  potential 
microcrack  will  be  less  than  a.  ctq  and  m  are  called  the  scale  and  shape  factors  respectively. 
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This  distribution  has  widely  been  used  in  analyzing  breakdown  phenomena^®’^®.  In  the 
simulations,  we  have  used  ctq  =  1.0  (thereby  mceisuring  all  stresses  in  units  of  cq)  and 
m  =  5.  One  observes  qualitatively  different  pictures  in  the  three  cases  shown.  In  the 
picture  corresponding  to  model  A,  we  find  a  highly  branched  macrocrack  (the  cluster  of 
cracks  that  connects  the  two  opposite  sides  of  the  lattice).  The  macrocrack  is  surrounded 
by  a  cloud  of  microcracks,  whose  width  seems  to  increase  in  the  direction  of  the  propagation 
of  the  macrocrack.  The  model  B  macrocrack  is  less  branched  and  the  microcrack  cloud 
around  it  seems  to  be  narrower  and  of  constant  thickness.  The  macrocrack  resulting  from 
model  C  simulation  is  topologically  almost  linear,  with  practically  no  branching.  There 
are  very  few  isolated  microcracks  around  the  main  crack. 

We  show  in  Fig.  2  the  nominal  stress  intensity  as  a  function  of  the  crack  length  for 
the  three  simulations  pictured  in  Fig.  1.  The  stress  intensity,  if/,  is  measured  in  units  of 
aoy/a^  a  being  the  original  crack  length.  (We  assume  that  during  the  whole  simulation  the 
condition  Aa  <  a  prevails.)  The  curve  for  the  model  A  simulation  exhibits  a  monotonic 
increase,  with  no  tendency  to  stabilize.  In  the  other  two  models,  the  nominal  stress 
intensity  needed  to  propagate  the  cracks  seems  to  saturate  at  a  stress  intensity  of  order 
unity.  The  final  stress  intensity  needed  in  the  model  A  simulation  is  about  1.7  (in  units  of 
ao\/a),  considerably  larger  than  this  asymptotic  stress  intensity.  The  different  behavior  is 
rooted  in  the  make  of  the  models:  in  model  A  each  stress  level  operates  only  once  while  in 
models  B  and  C  the  same  stress  keeps  forming  new  microcracks.  In  all  three  models,  there 
is  practically  no  crack  growth  below  a  stress  intensity  of  about  0.6,  since  the  probability 
of  having  a  potential  microcrack  whose  critical  stress  is  below  0.6ao  (  cf.  Eq.  (2))  is  very 
small. 

We  show  in  Fig.  3  the  total  number  of  microcracks  formed  as  a  function  of  the  increase 
in  the  length  of  the  macrocrack.  Again,  the  result  for  model  A  is  qualitatively  different  from 
the  other  two  models.  While  for  models  B  and  C,  N  increases  roughly  linearly  (constant 
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fracture  resistance),  the  increase  in  model  A  seems  to  be  with  a  power  larger  than  1, 
corresponding  to  an  increasing  R-curve.  We  shall  make  a  more  quantitative  comparison 
of  this  feature  later  on. 

Since  the  results  shown  in  Figs,  2  and  3  seem  to  have  some  noise,  reflecting  the 
stochastic  nature  of  the  lattice  (random  distribution  of  critical  stresses),  it  is  appealing 
to  simulate  a  large  number  of  lattices,  each  corresponding  to  a  different  realization  of  the 
distribution  of  the  critical  stresses,  in  order  to  obtain  smooth  results  which  Ccin  be  better 
analyzed.  We  simulated  200  different  realizations  of  the  three  models  described  above  with 
the  same  parameters  of  the  Weibull  distribution.  For  each  value  of  the  crack  length,  we 
averaged  the  corresponding  values  of  N  and  nominal  stress  intensity,  Kj.  Fig.  4  shows  the 
average  if/  as  a  function  of  Aa  of  the  main  crack.  The  stress  intensity  needed  to  sustain 
crack  growth  is  found  to  grow  rather  significantly  for  model  A,  but  is  almost  saturated  for 
models  B  and  C.  In  Fig.  5  we  show  N  as  &  function  of  the  crack  length  increment.  The 
inserts  in  Fig.  5  show  the  same  dependence  on  a  double  logarithmic  plot.  For  model  A, 
the  average  slope  in  the  range  5  <  Aa  <  60  is  about  1.4,  while  for  the  other  two  models, 
the  slope  is  about  1.  This  corresponds  to  a  flat  R-curve  for  models  B  and  C  and  to  an 
R-curve  that  is  approximated  as  Roc.  Aa®'^  for  model  A. 

So  far  the  results  have  been  for  the  case  in  which  the  screening,  as  defined  in  Eq.  (3), 
is  ignored.  Would  they  be  qualitatively  different  for  non-zero  values  of  the  screening 
pcirameter  a?  In  order  to  answer  this  question,  we  again  simulated  the  three  models  with 
the  same  Weibull  distribution,  using  a  =  0.6.  For  model  A,  the  same  qualitative  behavior 
is  found,  with  slightly  more  microcracks  around  the  main  crack.  Also  a  slightly  higher 
stress  intensity  was  needed  to  cause  the  same  crack  elongation.  This  indicates  additional 
toughening  due  to  microcrack  shielding,  although  the  same  R  oc  Aa®'^  is  maintained.  For 
model  B,  however,  a  qualitative  change  of  behavior  is  observed.  Fig.  6a  shows  iV  as  a 
function  of  Aa.  Using  a  double  logarithmic  plot,  we  find  that  the  behavior  is  consistent 
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with  N  oc  Therefore,  introduction  of  microcrack  shielding  changes  the  R-curve 

behavior  of  the  material  from  a  crack-length  independent  fracture  resistance  to  one  that 
increases  with  crack  elongation.  Furthermore,  when  a  is  increased  to  0.9,  we  find  an  even 
steeper  increase  of  the  resistance  to  fracture;  for  this  value  of  a,  the  exponent  increases 

from  1.2  to  about  1.4.  For  model  C,  the  situation  is  qualitatively  similar  to  the  case  with  no 

* 

microcrack  shielding.  Fig.  6b  shows  AT  as  a  function  of  Aa,  while  in  the  insert  we  present 
the  local  derivative  of  N.  We  see  that  the  derivative  is  a  slightly  decreasing  function  of 
Aa.  Note  that  this  decrease  should  stop  as  di\r/d(Aa)  reaches  the  value  of  1  as  at  least 
one  microcrack  is  needed  to  propagate  the  crack  by  one  lattice  unit. 

Another  interesting  question  is  how  the  R-curve  behavior  is  affected  by  the  distribution 
of  the  critical  stresses.  The  results  presented  so  far  have  been  for  a  Weibull  distribution 
of  critical  stresses.  We  have  repeated  the  simulations  for  models  B  and  C  for  different 
distributions.  In  the  extreme  case  in  which  all  the  potential  microcracks  have  the  same 
ffci  we  obtain  in  both  models  the  result  that  the  fracture  is  purely  brittle;  that  is,  the 
only  microcrack  formed  is  the  one  immediately  ahead  of  the  macrocrack  cind  the  creick 
propagates  at  the  critical  stress  level.  (In  this  case  there  is  no  randomness  in  the  system 
and  no  averaging  is  needed.)  This  result  is  not  surprising  since  the  most  stressed  bond  is 
always  the  one  immediately  in  front  of  the  crack.  We  also  tried  a  uniform  distribution  in 
the  range  ac,min  <  <^c  <  We  used  -  0.3  and  Cc^rnax  =  (^c,min  + 

case  the  results  are  similar  to  the  ones  obtained  with  Weibull  distribution;  an  increasing 
R-curve  behavior  is  observed  in  model  B  only  when  microcrack  shielding  is  significant. 
Otherwise  a  constant  resistance  to  fracture  as  a  function  of  crack  length  is  found.  For 
model  C,  the  resistance  to  fracture  remained  constant  (or  even  decreased  slightly  for  the 
quantity  defined  through  dNfd{Aa])  with  or  without  the  microcrack  shielding  effect  taken 
into  account.  For  this  model  we  also  tried  reducing  This  presumably  unrealistic 

limit  allows  for  the  formation  of  microcracks  at  arbitrarily  large  distances  from  the  crack 
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tip.  Probability  distributions  that  have  a  non-zero  value  for  vanishing  strengths  (or  the 
analogous  distribution  variable  in  a  different  problem)  are  known  to  lead  to  significantly 
different  critical  behavior^^*^®.  Even  for  this  limiting  distribution,  we  find  no  increasing 
R-curve  behavior  with  model  C,  although  both  the  number  of  microcracks  formed  and  the 
nominal  stress  intensity  at  any  crack  length  increment  are  larger  than  those  observed  with 
other  distributions  of  the  critical  stresses. 

IV.  DISCUSSION 

An  important  result  of  our  simulation  is  the  sensitivity  of  the  R-curve  behavior  to 
the  rate  of  loading  relative  to  the  internal  relaxation  processes  taking  place.  When  the 
loading  r.ate  is  higher  than  the  rate  of  stress  redistibution  in  the  material  or  the  inverse 
time  of  formation  of  a  single  microcrack,  we  find  an  R-curve  which  exhibits  a  significant 
increase  of  the  resistance  with  increasing  crack  length  (Figs.  4a  and  5a).  The  typical 
values  of  the  work  of  fracture  (or  of  the  nominal  stress  intensity)  observed  in  this  case 
(Model  A)  are  appreciably  higher  than  those  found  in  other  models  corresponding  to 
slower  rates  of  loading.  Such  high  rates  of  external  loading  are  hard  to  achieve  unless 
viscoelastic  effects  that  slow  down  the  internal  relaxations  are  present.  Indeed,  the  work 
of  fracture  in  a  typical  polymeric  material,  PMMA,  which  exhibits  viscoelastic  behavior, 
is  found  to  increase  linearly  with  the  crack  4relocity“®.  However,  the  fracture  mechanism 
in  this  case  is  through  craze  growth  governed  by  secondary  molecular  processes  and  not 
by  microcracking,  so  that  the  present  model  is  not  appropriate  for  such  materials.  To  our 
knowledge,  the  strain  rate  dependence  of  the  R-curve  in  a  material  with  an  established 
microcrack  toughening  mechanism  has  not  been  investigated  in  the  high  strain  rate  limit. 
As  our  model  indicates  an  enhanced  resistance  to  fracture  at  high  rates  of  loading,  it  would 
be  interesting  to  perform  an  experiment  to  check  this  prediction. 

The  other  two  models  that  we  have  defined,  models  B  and  C,  which  correspond  to  the 
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czse  in  the  extciml  Issuing  Is  s!zfxrej  Stan  ite  szlejczJ  pjtscesse^.,  fio  be  ii237£ 
T^srznt  to  ssksccjzcMzg  irgJtgsizSsj  which  ase  prei?ojsnai:j^  cSaslitizibei  than,  viscoehss- 
flc.  Tbs  sizin  mSssence  hetwreen  the  txo  mod^  Is  the  reSathre  izles  of  fl:e  tsio  isSeizal 
leSzxztlcs  pzocesses:  nasneSy,  the  relsxzilsn  of  tbe  Iscal  sixessas  axtd  tlie  foiH>2Sicxi  of  a 
smg!e  lolcrocTzck.  Model  B  is  appiopiizSc  wfeea  trse  foxxxiia&lon  of  gKoocracbs  izbes  place 
on  a  sborSer  tiaie  scale  than  tbe  local  stress  relaxa&ioa.  It  is  mieiesliiig  to  EoSe  that  tbe 
mode!  suggested  by  Hoaglaad  d  Eaibaiy*®  also  correspoads  to  the  time  r^me  inmlidt 
in  model  B,  altboa^  no  mention  of  time  scales  is  made  in  tbeir  work.  The  reason  for 
this  obserraiion  is  that  in  their  model,  as  in  onr  model  B,  all  tential  microcrack  sites  on 
which  the  (self-consistently  determined)  stress  exceeds  the  critical  stress  are  allowed  to  be¬ 
come  microcracks.  This  implicitly’  assumes  that  the  local  stresses  do  not  react  fast  enou^ 
to  the  formation  of  microcracks,  since  othenvisc  the  formation  of  a  single  micio,rack  (the 
one  with  the  largest  excess  stress,  as  in  our  model  C)  would  change  the  stress  field  in  the 
sample. 

We  should  note  that  of  the  three  models  only  model  C  is  really  quasi-static  since  the 
order  in  which  microcracks  are  formed  is  independent  of  the  increment  of  the  external 
stress.  In  the  simulations  we  performed  we  increased  the  nominal  stress  intensity  in  steps 
of  0.05.  For  some  cases  we  repeated  the  simulations  with  half  the  step  size  and  obtained 
identical  R-curves  for  model  C.  As  one  can  expect,  model  A  was  especially  sensitive  to 
this  step  size.  Both  measxires  of  the  crack  resistance  exhibited  a  decrease  with  decreasing 
step  size;  the  qualitative  behavior,  however,  did  not  change.  The  nominal  stress  needed  to 
sustain  crack  growth  increased  as  previously.  The  number  of  microcracks  formed  at  creick 
length  Ac  grew  also  proportionally  to  Aa®  as  with  the  larger  step  size.  Although  not 
inherently  quasi-static,  model  B  exhibited  no  detectable  difference  between  the  R-curves 
obtained  with  both  step  sizes.  This  shows  that  by  allowing  the  stresses  to  relax  after  the 
formation  of  microcracks  and  before  increasing  the  nominal  stress  intensity,  we  compensate 
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fcx  ihe  fgnssa^-sa  m  J2®ie  ihzn  cze  iSLCsocszck  sataSaneccs^  zszd  ike  exhSb'Hs  a 
gszsi-staSic  b^sawg. 

Hesalis  for  03r  B  sbov  feat  aa  iscreasisg  R-cnnre  is  observed  oidf  wbea  fee 

critical  stresses  of  fee  potea&I  rnicrocracks  has  some  raadom  distiibiifem  (^gs.  4b  and 
5b).  Whea  all  fee  critical  stresses  hasre  fee  same  valae  with  no  dispersion,  we  obtain  a 
constant  hzdtaie  energy,  ^ea  by  fee  energy  of  formation  of  a  single  microcrack.  If  we 
nse  fee  alternative  deSnition  of  fee  frzctme  resistance  whlfe.  is  related  to  the  nominal 
stress  intensity,  this  case  corresponds  to  an  unstable  crack  propogation  at  a  constant 
stress  intensity,  after  a  transient  rise.  When  the  critical  stresses  have  some  distribution, 
then  fee  results  depend  upon  the  intensity  of  interaction  between  the  macrocrack  and  fee 
microcracks.  If  we  ignore  the  crack-microcrack  interactions,  the  number  of  microcrac^ 
formed  to  increase  the  crack  length  by  an  amount  Aa  is  found  to  grow  linearly  •  Uh  Ao, 
indicating  a  constant  resistance  to  fracture.  The  nominal  stress  intensify  to  sustain  that 
crack  propagation,  on  the  other  hand,  shov/s  a  mild  increase  with  crack  length.  This 
discrepancy  between  the  two  dciinitions  of  the  resistance  to  fracture  has  already  been 
found^®  in  ceramic  materials  undergoing  what  b  termed  ‘transformation  toughening’^. 
This  is  a  toughening  mechanism  due  to  the  transformation  of  zirconia  particles  dispersed 
in  some  ceramic  materials  like  AI2O3  or  Si3N4  from  a  tetragonal  symmetry  to  monoclinic 
synunetry  under  the  stress  field  of  the  m2u:rocrack.  The  transformation  process  acts  as  an 
additional  dissipative  energy  sink,  thus  increasing  the  energy  input  needed  to  propagate 
the  main  crack.  Rose  and  Swain^®  found  that  the  two  definitions  of  the  resistance  show  a 
qualitatively  different  transient  behavior-  the  quantity  defined  through  energy  decreases  to 
a  steady  state  value,  while  the  stress  intensity  based  definition  increases  to  the  same  value. 
In  our  simulations,  the  discrepancy  is  less  dramatic  but  still  noticeble.  When  the  crack- 
microcrack  interaction  is  taken  into  account,  the  number  of  microcracks  formed  increases 
more  rapidly  than  the  crack  length  increment,  giving  an  increasing  resistance  to  fracture. 
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The  vzctezse  in  the  resssience  ts  fonnd  to  he  ht^er  the  stronger  mxcrocxack  ihe  slueSding 
deSsed  in  Eo.  (3). 

Li  tMe  xe^me  of  time  scales  for  whidi  model  C  is  appropriate,  we  find  that  microc- 
xacking  does  not  resnlt  in  a  rising  R-curve  bdiavior,  even  when  microcrach  shielding  is 
tahen  into  acconnt  and  a  wide  distribution  of  critical  stresses  are  present  in  the  materia! 
(although  the  stress  intensity  as  a  function  of  the  crack  length  esdiibits  a  slight  increase, 
as  in  model  B.).  This  is  a  significant  result  as  w'e  think  that  model  C  is  relevant  to  some 
realistic  cases,  since  we  expect  the  microcrack  formation  time  to  he  comparable  to  or  larger 
than  the  relaxation  time  of  the  local  stress  distribution  around  the  macrocrack  for  elastic 
materials.  To  obtain  a  qualitiative  estimate  for  these  two  time  scales,  let  us  note  that 
in  an  elastic  material,  the  local  stress  relaxation  will  propagate  with  the  material  sound 
velocity,  which  is  of  order  10^msec“^.  The  stresses  in  a  region  of  order,  say,  1cm -will 
relax  within  a  time  of  order  r/  10"®  sec.  Using  acoustic  emission  data  on  fracture  in 
concrete^'^^,  the  duration  of  the  microcrack  formation  process  can  be  estimated  to  be 
10~®  sec  or  longer.  This  is  the  same  order  of  magnitude  as  r;  and  model  C  may  be  relevant 
in  some  materials.  Furthermore,  several  models  proposed  for  some  breakdown  phenonema 
(electrical^®*^^  as  well  as  mechanical^^*^^)  also  assume  this  time  regime.  Our  results  in¬ 
dicate  that  in  such  cases,  the  resistance  to  fracture  will  be  constant  as  i  function  of  the 
crack  length  increment. 
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V.  CONCLUSIONS 

We  hzve  described  3iree  simple  models  of  crack  propagation  through  microcrack  coales¬ 
cence  under  uniaxial  tmirile  loading.  The  models  difier  in  the  relative  time  scales  of  the 
plQrsIcal  processes  that  occnr  in  the  ^stem.  We  find  that  when  the  rate  of  loading  is 
faster  than  the  iniemal  relaxations  in  the  ^'stem  (microcrack  formation  and  local  stress 
relaxation  following  any  change  in  the  ^stem),  the  microcracking  material  exhibits  an 
increasing  resistance  to  fracture  with  increasing  crack  length.  We  do  not,  however,  expect 
this  model  to  be  relevant  for  microcracking  elastic  matcriak  like  ceramics  since  it  is  diScult 
to  realize  such  high  loading  rates.  When  the  time  scale  involved  in  local  stress  relaxation 
is  shorter  than  the  inverse  rate  of  loading,  but  longer  than  the  time  of  formation  of  a  single 
nucrocrack,  increasing  fracture  resistance  behavior  is  observed  only  if  the  material  has  a 
distribution  of  critical  stresses  of  potential  microcracks  and  if  there  is  substantial  microc¬ 
rack  shielding.  When  the  interaction  between  the  crack  and  the  microcracks  are  such  that 
they  tend  to  cancel^^,  a  brittle  behavior  with  a  constant  resistance  to  fracture  is  observed. 
When  the  local  stress  relaxation  is  the  fastest  process,  we  find  that  microcracking  does  not 
lead  to  any  toughening;  that  is,  the  resistance  to  fracture  does  not  increase  with  increasing 
crack  length. 

In  summary,  this  work  has  shown  the  importance  of  microcrack  shielding  to  a  rising  R- 
curve.  More  and  better  data  on  shielding  (both  experimental  and  theoretical)  will  justify 
further  computations  along  the  lines  proposed  in  this  study. 
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FIGURE  CAPTIONS 

Figure  1.  The  macrocrack  and  the  array  of  microcracks  formed  in  the  system  when 
the  macrocrack  has  traversed  the  lattice  of  linear  size  L  =  60.  (a)  Model  A,  (b)  Model  B 
and  (c)  Model  C.  In  all  the  cases  a  Weibull  distribution  of  critical  stresses  with  ojo  =  1.0 
and  m  =  5  was  used.  Crack-microcrack  interaction  parameter  is  set  to  0. 

Figure  2.  The  nominal  stress  intensity,  if/,  needed  to  increase  the  macrocrack  length 
by  Aa'for  the  three  structures  shown  in  Fig.  1. 

Figure  3.  The  total  number  of  microcracks  formed  (including  those  that  have  joined 
the  microcrack)  in  the  three  structures  shown  in  Fig.  1,  when  the  macrocrack  length 
increased  by  A  a. 

Figure  4.  The  average  if/  as  a  function  of  crack  length  increment  for  (a)  Model  A,  (b) 
Model  B  and  (c)  Model  C.  The  average  is  over  200  realizations  of  a  60  X  60  lattice  with 
a  Weibull  distribution  of  critical  stresses,  with  oq  =  1.0  and  m  =  6.  Crack-microcrack 
interaction  parameter  is  set  to  0. 

Figure  5.  The  average  number  of  microcracks  formed  as  a  function  of  crack  length 
increment  for  (a)  Model  A,  (b)  Model  B  and  (c)  Model  C.  The  average  is  over  200  realizar 
tions  of  a  60  X  60  lattice  with  a  Weibull  distribution  of  critical  stresses,  with  (Tq  =  1.0  and 
m  =  5.  Crack-microcrack  interaction  parameter  is  set  to  0.  The  inserts  show  the  same 
dependence  on  a  double  logarithmic  plot. 

Figure  6.  The  average  number  of  microcracks  formed  as  a  function  of  crack  length 
increment  for  (a)  Model  B  and  (b)  Model  C.  The  average  is  over  200  realizations  of  a  60x60 
lattice  with  a  Weibull  distribution  of  critical  stresses,  with  oq  =  1.0  and  m  =  5.  Crack- 
microcrack  interaction  parameter  is  set  to  0.6.  The  insert  in  (b)  is  the  local  derivative. 
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Abstract.  After  introducing  the  formalism  of  maximum  entropy  and  reviewing  alternative 
approaches  for  fragment  size  statistics,  this  paper  derives  a  general  distribution  law  (similar 
in  form  to  the  Bose-Einstein  statistics)  and  applies  it  to  distributions  observed  in  rock  mining, 
exploding  metallic  shells,  shattered  crystal  pieces,  droplets  in  spray,  atomic  or  molecular 
clusters,  space  debris  and  fragmented  nuclei.  Variations  of  fragment  number  with  size 
that  are  power-law-like  (fractal),  humped  or  exponential  can  lead  to  physically  significant 
conclusions  regarding  the  fracturing  mechanism.  Theoretical  aspects  of  the  maximum- 
entropy  method  in  the  derivation  of  the  distribution  law  (including  some  inherent  difficulties) 
are  discussed. 


1.  Introduction 

The  pieces  of  a  broken  object  are  among  the  most  common  examples  of  randomness 
and  disorder  (Zurek  1989 ,  section  V) .  There  may  be  psychological  reasons  for  associating 
disorder  with  broken  pieces,  perhaps  due  to  our  dismay  at  seeing  a  nicely  designed 
object  shattered  and  to  the  need  to  clean  up  the  resulting  mess,  but  it  is  a  fact  that  it  has 
been  found  difficult  to  read  some  regularity  into  fragments.  The  regularity  that  we  can 
hope  to  get  is  a  law  describing  the  distribution  of  fragments,  especially  with  respect 
to  their  sizes  but  also  regarding  shapes  and  internal  parameters.  The  prediction  of 
distribution  has  in  many  cases  practical  advantages,  relevant  to  scientific  and  tech¬ 
nological  issues.  A  few  instances  taken  from  physics  will  be  given  later  (section  6);  at 
this  point  I  shall  describe  the  industrial  importance  of  oil-shale  fragment  size  distribution 
since  it  was  through  this  that  I  became  involved  in  fragmentation. 

In  a  standard  mode  of  operation  (called  ‘bench  blasting’)  the  bituminous  rock  (oil 
shale)  is  explosively  mined  to  produce  the  debris  containing  fragments.  The  large 
specimens  are  then  comminuted  mechanically  to  bring  them  down  to  a  size  such  that 
retorting  (or  extraction  of  oil  by  partial  burning)  can  be  carried  out.  Since  comminution 
takes  time  and  money,  large  fragments  (exceeding  about  50  cm  in  diameter)  are  unde¬ 
sirable,  and  so  are  small  pieces  (less  than  0.5  cm)  since  these  clog  up  the  retort.  The 
mining  has  to  be  planned  accordingly  and,  while  one  rarely  has  control  over  the  geo¬ 
mechanical  properties  (faults,  cracks,  etc)  of  the  rock,  one  can  still  vary  some  parameters 
of  the  mining,  namely,  location  of  explosion  boreholes,  temporal  sequence  of 
their  activation  and  type  of  explosive.  By  correct  choice  of  these  one  can  hope  to  obtain 
a  fragment  size  distribution  with  the  desired  properties  of  not  having  too  many  pieces 
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outsde  the  lower  and  upper  Emils.  The  task  ai  hand  is  ths2s  the  derivsiiiKi  of  a  feagment 
size  disfribuEon  in  the  context  of  a  phy^cal  process  {the  passage  of  detcnaiion  sliess 
waves)andinagh'engecinetricseitiEg(someborefeo!ef®iferaajEdi3aa]orfes2ltsssiiitose 
in  the  rock)  (Jaeger  et  al  1986b)-  At  the  same  time,  several  physcal  factors  will  remain 
unknown,  like  subsurface  joints  and  cracks,  leaving  Eie  theoretician  with  an  under- 
defined  assignmenL 

It  w'ould  appear  that  the  mrcdmum-entropy  method  (dSESi)  is  an  ideal  tool  to  hand!e 
this  situation  and  others  in  which  we  wish  to  predict  fragment  distrfmiiioiis.  Reasons 
both  favouring  and  opposing  this  eiq^ctation  are  outlined  in  the  next  sectk)ii  and  in 
appends  1,  where  we  return  to  a  discussion  on  the  clmms  of  iJESk!. 

Fracture  of  soEds  as  a  ph}^ca]  process  is  the  subject  of  appendix  2. 


2.  The  method  of  entropy  maximization 

Let  us  denote  by /the  events  that  we  wish  to  predict  (like  the  event  of  finding/*  feagments 
of  a  ^ven  size  in  a  debris  heap)  and  by  p/the  corresponding  probabiEties.  We  then  assert 
that  the  mem  rives  the  most  conservative  probabiE^  distribution  p,  subject  to  the 
available  knowledge,  which  enters  the  formalism  as  ‘constraints’.  The  maximum- 
entropy  solution  is  that  which  maximizes  the  information  entropy  (also  called  ‘misring 
information’) 

5=-Sp/logp/  (1) 

/ 

(Shannon  1948)  subject  to  the  constraints.  The  constraints  are  assured  mathematically 
by  Lagrange  multipliers  (appendix  1). 

The  meaning  of  ‘most  conservative’  in  the  assertion  of  the  last  paragraph  is  twofold. 
First,  the  mem  probabilities  Pj  are  about  the  smoothest,  flattest  functions  of  the  events/ 
(suitably  ordered)  consistent  with  the  constraints;  secondly,  the  predictions Pj  make  do 
with  just  the  constraints  explicitly  imposed  and  not  with  any  other  constraints  perhaps 
tacitly  implied  (in  obtaining  non-MEMp^).  The  latter  meaningprovides  asimple,  heuristic 
proof  of  the  assertion.  It  can  be  shown  that  imposing  any  additional  constraints  reduces 
S,  i.e.  it  reduces  the  ‘missing  information’.  If  our  choice  for  the  probabilities  Pj  is  non- 
MEM,  this  choice  will  reduce  S  below  its  maximum  value  and  will  imply  some  extra 
constraint.  Therefore  the  mem  solution  is  the  only  one  consistent  with  the  explicitly 
stated  knowledge. 

More  impressive  claims  have  also  been  made  for  the  mem  as  a  predictive  method 
(Jaynes  1983,  Skilling  1984,  Rosenblatt-Roth  1987,  1988),  In  the  context  of  frag¬ 
mentation,  one  can  show  that,  if  N  independent  fragmentation  experiments  are 
performed,  the  mean  of  the  size  distributions  will  tend,  asymptotically  as  to  the 

MEM  value  of  Pj. 

A  consistency  criterion  was  enunciated  by  Tikochinsky  et  al  (1984)  in  the  following 
sense:  the  distribution  of  mem  is  the  only  one  that  will  be  confirmed  by  successive 
repetition  of  experiments. 

In  view  of  the  existence  of  several  time-honoured  and  practised  predictive  procedures 
(like  maximum-likelihood,  least-squares  fitting,  Bayesian  predictions),  the  preceding 
claim  is  bound  to  have  raised  opposition  (Titterington  1984).  To  answer  this.  Skilling 
(1984)  writes;  ‘I  find  the  arguments  for  using  maximum  entropy  to  be  deeply  compelling. 
[mem]  stands  in  splendid  isolation,  on  fundamental  grounds.’ 
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But  evea  tiiis  is  not  all,  scce  acooidii^to  Jaynes  (19S3),  mem  predicdoas  cannoihe 
lEssiakea:  if  ihsy  are  cot  sopported  by  expenment,  they  aiust  be  regarded  as  disclo^g 
(the  need  for  or  ensfence  o^  addiliana]  infonnational  constraints.  Jaynes  is  to  be 
credited  vsith  the  use  of  iofonnation  entropy  as  a  predictive  tool  in  phjscs  and  this 
review  does  perhaps  no  more  than  describe  the  implementation  of  this  pro^amme  for 
fragmentation.  In  his  works  (collected  in  Jaynes  (19S3))  thermodynamic  entropy  is 
represented  as  a  measure  of  the  insuSdent^*  in  our  knowledge  of  a  phydcal  ^■stem.  In 
an  op|X)site  view  (Denbigh  and  Denbigh  198^  thermodynamic  entropy  is  regarded  as 
an  objective  attribute  of  the  system.  In  a  compromise  approach  (Zurek  1989)  the 
thermodynamic  entropy  is  a  sum  of  (X)ntributions  of  phydcal  randomness  and  of  the 
observer’ s  ignorance. 

Without  taking  ades  in  this  controversy,  we  wish  to  demonstrate  here  that  the 
maximum-entropy  soludon py  with  energy  E ^'en  is  formally  identical  to  the  minimum- 
energy  solution  Pj  with  entropy  S  ^'en.  With  the  introduction  of  Lagrange  multipliers, 
denoted  by /5,  V  and  r,u(respectively)forthetwoproblems,  the  quantities  (henceforth 
called  Lagranaans)  to  be  extremized  are 

L  =  S-p{E-E)  =  -  ^pj  logpj  -  ^^2  EjPf  -E^-  v^pi  -  (2) 

and 

-  A  =  £  -  T(S  -  S)  =  2  £/P/  +  Pi  log  Py■^  -  #1(2  Pi  - 1).  (3) 

The  sets  of  equations  dLjdpj  =  0  and  dAfdPj  =  0  have  evidently  formally  the  same 
solutions.  These  become  equ^. 

Pi  =  =  exp(-^/r)  j  2  ExpirEilT)  (4) 

if  one  supposes  the  relations 

r=^-i  =  £/5.  (5) 

Clearly  equation  (4)  is  also  the  solution  obtained  by  minimizing  the  free-energy 
expression  (for  temperature  T  given).  The  maximum-entropy  procedure  would  thus 
appear  to  be  beyond  doubt,  but  questions  remain  whether  the  system  (especially  a 
classical  one)  is  truly  ergodic  over  the  events  /. 


3.  Distribution  of  fragments 

3.1.  Exponential  distributions  based  on  mem 

Apparently  the  earliest  use  of  mem  for  fragment  size  distribution  is  due  to  Griffith  (1943) . 
He  supposed  that: 

(i)  a  given  quantity  of  energy  is  converted  to  fragmentation  (this  quantity  depending 
on  the  nature  and  intensity  of  the  fracturing  process); 

(ii)  the  energy  goes  entirely  into  the  creation  of  free  surfaces;  and 

(iii)  each  molecule  in  the  solid  competes  for  the  available  energy  on  equal  footing. 
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If  Ng  denotes  the  number  of  molecules  belon^g  to  5-sized  ^agments.  Njs  is  the 
number  of  s-sized  fragments.  The  energy  invested  in  each  fragment  is  where  y  is 

the  surface  energy  density  and  g  a  geometrical  factor.  The  constraint  is 

£  =  yg  2  =  yg  2  (6) 

5  S 

so  that  identifying  Sj  in  equation  (2)  with  Ygs~'^  and  regarding  the  solutions  in  equation 
(4)  one  derives 

oc  expC-A/s’^)  (7) 

where  A  is  a  constant  to  be  evaluated  from  equation  (6).  In  practice,  in  the  break-up  of 
solids,  the  number  of  fragments  Njs  decreases  with  size  much  faster  than  that  ^ven  by 
equation  (7),  indicating  that  assumption  (iii)  places  too  much  weight  on  large  fragments 
and  that  sinks  of  energy  other  than  the  surface  energy  (involved  in  assumption  (ii))  are 
also  present.  Moreover,  molecules  in  the  middle  of  the  fragment  do  not  feel  the  fracture 
process  in  the  same  way  as  do  molecules  on  the  surface .  We  shall  see  that  when  fragments 
are  formed  by  accretion  of  molecules,  assumption  (iii)  has  an  empirical  basis  (section 
6.4). 

In  a  discussion  of  the  size  distributions  of  exploding  shells  and  of  other  two-dimen¬ 
sional  objects,  Grady  and  Kipp  (1985)  used  only  a  geometrical  constraint,  namely 

A=g'2n^s  (8) 

S 

meaning  that  the  number  of  plane  fragments  of  size  s  times  their  area  gs  (where  g  is 
again  a  geometrical  factor)  sums  up  to  the  total  area  A.  The  MEM  soution  is  clearly 

n,  oc  exp{-gsjA).  (9) 

This  is  to  be  contrasted  with  the  semiempirical  distribution  due  to  Mott  (1947)  (see  also 
Mott  and  Linfoot  1943)  in  which  the  cumulative  distribution  arising  from  is  an 
exponential  function  of  the  linear  size,  namely  in  two  dimensions  and  in  three 
dimensions.  Data  on  metallic  shells  have  been  accounted  for  by  Mott’s  distribution 
(Sternberg  1973),  though  discrepancies  appear,  which  are  treated  in  greater  detail  in 
section  6.1.  Deviations  from  Mott’s  law  have  been  found  in  fragmented  rock  data  and 
in  a  numerical  simulation  (Englman  et  al  1984). 

5.2.  Fragmentation  by  a  fast,  sudden  process 

Fracture  by  fast  energy  input,  such  as  occurs  in  explosions,  projectile  or  laser  impact,  a 
shattering  event,  explosive  expansion  or  a  sharp  temperature  change,  has  been  the 
subject  of  dynamic  theories  (Mott  1947,  Grady  1982,  Grady  et  al  1985,  Glenn  and 
Chudnovsky  1986,  Glenn  etal  1986,  Curran  etal  1987). 

These  theories  are  incorporated  in  the  probabilistic  prediction  of  fragmentation  that 
we  describe  now.  Characteristic  of  fast  processes  is  that  fracture  events  take  place  at 
different  points  of  the  solid  independent  of  one  another.  This  is  because  newly  created 
fractures  have  no  time  to  communicate  with  each  other ,  except  that  the  energy  invested  in 
the  fracture  cannot  exceed  the  total  energy  available  for  fracture  and  that  the  fragmented 
volumes  are  also  bounded  by  the  original  volume.  Following  Grady  (1982),  we  suppose 
that  all  discontinuities  formed  by  the  stress  wave  (e.g.  cracks,  voids)  immediately 
become  faces  of  fragments.  An  alternative  view,  which  remains  to  be  explored,  would 
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be  to  suppose  that  in  the  firet  place  the  fast  process  creates  fractures  that  can  then  join 
up  to  enclose  a  fragment.  The  advantage  of  this  second  view  is  that  unfinished  fractures, 
called  ‘internal  damage’  (Jaeger  et  d  1986a,  c,  Yatom  and  Ruppin  1989)  are  also 
envisaged  to  be  present  as  the  result  of  stress. 

Following  a  more  detailed  account  in  the  original  paper  (Englman  et  d  1988b)  we 
associate  with  each  fragment  of  linear  size  a  an  energy  e(a)  (see  equation  (10)  below) 
which  in  essence  contains  all  the  physics  entering  the  formation  of  fragments  (figure  1). 
Summing  e(a)  over  all  fragments,  we  equate  the  sum  to  the  total  energy  input  E  for  the 
fragmentation  process;  this  provides  the  firet  constraint  (equation  (15)  below) ,  Implicit 
in  this  procedure  is  the  assumption  that  at  the  instant  of  fragmentation  all  secondary 
featuresoftheprocess(suchasheat,noise,fragmentmotion)stillresideinthefragments. 
Our  way  of  invoking  energy  conservation  differs  from  that  of  Glenn  et  al  (1986),  who 
balance  the  created  surface  energy  against  the  kinetic  and  stored  energies.  The  present 
approach  appears  justified  (barring  the  role  of  expanding  gases  formed  in  an  explosion) 
because  it  includes  the  physical  state  of  the  fragments  immediately  after,  fragmentation. 
We  assume  fragments  of  linear  size  a  and  roughly  spherical  shape  whose  energy 

e(d)  =  (2n:/3)pu^a^  +  (2jcf5)pPa^  +  (4jtf6K)a~a^  +  4;ryfl^(l  -  a/R)  (10) 

consists  of  the  following  terms:  the  kinetic  energy  (ke)  of  the  fragment  of  mass  density 
p  moving  as  a  rigid  body  with  speed  u  (this  term  takes  up  the  largest  part  of  the  energy); 
the  internal  ke  about  the  centre  of  mass  (expressed  in  terms  of  e^,  the  volumetric 
dilatation  rate  squared,  or  the  trace  of  the  square  of  the  stretching  tensor);  the  stress 
energy  stored  in  the  fragments  (Glenn  et  al  1986)  including  sound  and  elastic  waves  (o 
being  the  stress  and  K  the  bulk  modulus);  and  the  surface  fomiation  energy  (y  is  the 
surface  energy  density  and  R  is  the  radius  of  the  spherical  solid  before  fragmentation). 
Other  terms  depending  on  size  a  or  on  other  variables  (e.g.  shape  or  chemical  compo¬ 
sition)  may  be  added  to  equation  (10). 

The  probability  of  having  fragments  of  size  a  is  written 

p(n„,a)  with  Sp(nfl,fl)  =  l  (11) 

tla 

and  the  most  likely  number  of  such  fragments  and  their  masses  are,  respectively, 

(«fl>  =  EnaP(nfl,o)  (12) 

tla 

and 

w(a)  =  (4n:/3)pa^(«a).  (13) 

Assuming  that  a  fixed  mass  W  of  materials  is  being  fragmented,  we  have  the  constraints 
W=Sw(fl)  =  W  (14) 

a 

and 

£  =  S  e(p){n.)  =  t  (15) 

a 

where  the  second  constraint  is  due  to  energy  conservation. 

(In  equation  (11),  we  have  chosen  as  the  fundamental  event  whose  probability p(no, 
a)  is  sought  the  finding  of  ria  fragments  of  linear  size  a,  rather  than  the  event  of  finding 
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at  random  an  a-sized  fragment,  whose  relative  frequency  is p{a).  Our  choice  reflects  the 
practice  in  mineral  processing  of  the  operation  of  sieving.) 

In  the  sense  of  the  MEM,  the  best  distribution  p  subject  to  the  constraints  is  obtained 
by  maximizing 

L  =  -  2  S pK,  a)  logp(«o,  a)-p{E-E)-  p{W  -  W)  (16) 

a  rta 

through  variation  of  p{na,  a)  and  adjusting  the  Lagrange  multipliers  P  and  p  so  that  the 
constraints  (13)  and  (14)  are  satisfied. 

The  solution  is 


pin^,a)  =  exp(-n^x„)/2  exp(-/2^xj 

tta 


(17) 


where 

Xfl  =  pe(a)  +  p{4jtf3)pa^  (18) 

and 

<n,>  =  (e^«-l)-‘  (19) 

(20) 

for  the  majority  of  fragments,  i.e.  excluding  large-sized  fragments,  for  which  >  1, 
whose  mean  number  decreases  exponentially.  In  (18)  the  fragment  mass  and  energy  are 
both  polynomials  in  a,  with  powers  0  ranging  from  2  to  5  (equation  (10)) .  In  a  size  range 
where  a  single  power  law  dominates  in  (18),  the  distribution  of  sizes  will  exhibit  a  power- 
law  behaviour  in  the  form 

I  da{na)°^a~^  (21) 

a 

with  the  ‘fractal  dimension’  D  being  given  in  terms  of  the  predominant  values  of  0  as 
D  =  0-1. 

Turcotte  (1986)  has  listed  21  objects  whose  fragment  size  distribution  behaves 
according  to  (21),  with  D  falling  between  1.5  and  3.5,  i.e.  conforming  to  the  theoretical 
values  in  equations  (10),  (18)  and  (20).  Experiments  on  simulation  of  astroid  impacts 
gave  D  ~  2.5  (Capaccioni  et  al  1986). 

To  determine  the  parameters  ^  and  p  in  the  distribution  (n)  (equations  (18)-(20)) 
from  the  constraints,  we  first  combine  (14)  and  (15)  to  write 

/  =  i8E-t-piy  =  2x,(e*--l)-*  (22) 


where  the  summation  is  over  all  sizes  that  are  collected  in  discrete  bins.  Let  the  number 
of  bins  (or  summands)  be  4*  then 

l<h  (23) 

since  the  summand  is  less  than  unity.  The  dominant  contribution  to  e{a)  and  Xg  is  a  term 
proportional  to  the  fragment  volume  a^.  Thus,  1  is  proportional  to  W. 
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Figure  1.  Forces  acting  on  a  fragment.  The  figure 
illustrates  the  body  forces  pu  (inertial)  and  a 
(uniform  stress),  the  surface  energy  density  y 
and  the  transverse  dynamic  strain  e. 


Figure  2.  Weight  distribution  of  fragments, 
equation  (24).  The  weight  iv(fl)  and  the  linear 
fragment  size  a  are  normalized  to  their  values  at 
a  (equation  (26)).  Beyond  the  dis¬ 

tribution  drops  off  exponentially  with  a  and 
equation  (24),  on  which  this  figure  is  based,  is 
not  applicable.  The  parameter  v  in  equation  (27) 
has  the  value  10“^. 


The  weight  distribution  (13)  for  the  majority  of  fragments  (equation  (10))  is 
expressed  in  a  reduced  form  as 

iv(n)  =  W- V(2w-‘  +  1  +  us^)  (24) 

s  =  afd.  (25) 

It  has  a  maximum  at  a  =  a,  given  by 

a  =  (5y/p£2)i/3  (26) 

which  serves  as  a  convenient  scaling  parameter  for  the  fragments  (Grady  1982).  The 
other  parameter  v  in  the  distribution  (24)  contains  the  multipliers  j8  and  p,  which  must 
be  eliminated  by  calculating  the  mass  (or  weight)  constraint  (14).  The  result  is 

V  =  {Tcf3){dlll)^  (27) 

where  I  is  the  bin  width,  so  that,  as  shown  in  equation  (22),  II  is  the  effective  size  range 
scanned  by  the  sampling  process.  The  saddle-point  d  must  fall  well  inside  the  size  range, 
so  that 


V  (28) 

The  weight  distribution,  shown  in  figure  2,  has  a  maximum  at  a  =  a,  is  broad,  and  is 
asymmetric.  Thus,  the  maximum  derived  by  mem  coincides  with  the  characteristic 
fragment  size  of  Grady  (1982)  and  depends  (through  equation  (26))  only  on  the  dynamic 
fracture  parameters  y  and  e^.  On  the  other  hand,  a  is  independent  of  either  the  centre- 
of-rnass  ke  or  of  the  stored  stress  a,  or  indeed  of  any  contribution  to  the  energy  that  is 
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proportional  to  the  fragment  volume  4;rcV3-  The  mathematical  explanation  for  this  is 
that  in  the  summand  of  (14), 

w(fl)  cc  a^ria  ~  a^x~^  =  [^e(a)fa^  +  ju]"* 

all  such  contributions  appear  as  an  additive  constant  to  the  Lagrange  multiplier  fi,  which 
is  a  parameter  to  be  adjusted. 

The  empirical  role  of  d  is  well  documented  in  the  fragmentation  of  oil  shale  (Grady 
1982),  brittle  steel  (Weiner  and  Rogers  1979)  and  rapidly  heated  water  (Blink  and 
Hoover  1985) .  Since  we  find  that,  in  general,  the  mean  size  averaged  over  the  distribution 
exceeds  its  maximum  d,  a  more  careful  determination  of  the  average  and  a  specification 
of  the  averaging  procedure  are  called  for. 

In  ‘cratering’  situations,  where  only  part  of  the  solid  is  fragmented,  the  constraint 
(14)  on  the  total  weight  W  does  not  apply  and  (i  does  not  enter.  Yet,  remarkably, 
elimination  of  the  parameter  leads  to  the  same  distribution  (24)  as  in  the  fixed  W 
situation.  This  was  subject  to  test  when  fused  alumina  cement  targets  were  partially 
(—30%)  and  almost  fully  (90%)  fragmented  by  fast  aluminium  projectiles  (Bianchi  et  al 
1984) .  The  fractal  exponents  were,  respectively,  D  =  2.76  and  3.06,  which  do  not  support 
our  result,  though  their  magnitude  is  within  the  range  predicted  in  equations  (10),  (18) 
and  (20). 

The  weight  W  of  the  crater  mass  scales  with  the  energy  input  E  according  to 

W  =  {2jtpEd/9Y)idlll)\  (29) 

The  threshold  energy  input  E,  below  which  fragmentation  cannot  occur  (the  con¬ 
straint  equation  (15)  having  no  solution)  is  given  by 

E^  =  Anyd^{d/l).  (30) 

A  different  threshold  energy  for  fragmentatiqn,  resulting  from  a  geometrical  constraint, 
was  obtained  by  Yatom  and  Ruppin  (1989). 

A  striking  result  of  the  mem  distribution  is  the  prediction  of  fractal  distribution, 
equations  (17)-(20),  with  fractal  exponents  in  the  range  of  the  observed  ones  (Turcotte 
1986).  The  result  depends  on  two  circumstances;  first,  on  the  property  that  the  energy 
and  volume  terms,  given  by  e{a)  and  go  to  zero  as  the  size  a  vanishes  so  that  there  is 
no  contribution  to  the  constraints  from  o-^O,  admitting  of  a  ‘Bose  condensation’; 
secondly,  on  the  absence  of  a  constraint  on  the  total  number  of  fragments.  Since 
the  latter  circumstance  is  tied  to  the  non-atomistic  description,  it  seems  that  fractal 
distributions  might  be  the  rule  for  classical  objects.  We  emphasize  though  that  fractal, 
inverse  power-law  behaviour  is  predicted  only  for  the  low  size  end;  at  large  sizes  an 
exponential  decrease  with  size  takes  place.  In  our  discussion  of  experimental  dis¬ 
tributions  in  the  sequel,  the  question  of  power-law  versus  exponential  form  will  come 
up  frequently. 

In  summary  of  this  section,  an  essentially  simple,  two-parameter  (a  and  v)  dis¬ 
tribution  (equation  (24))  has  been  derived,  based  on  the  broad  principle  of  maximum 
entropy  and  incorporating  the  physics  of  fragmentation  through  the  energy  of  fragments 
of  size  a  (equation  (10)).  The  typical  Grady  fragment  size  has  been  rigorously  rederived 
and  identified  with  the  saddle-point  d  and  the  maximum  of  the  distribution;  d  was  found 
to  be  unaffected  by-mechanisms  that  are  coextensive  with  the  fragment  volume. 

5.5.  Statistics  of  propagating  fragmentation 

In  the  previous  section  we  have  described  a  fragmentation  process  that  is  fast,  sudden 
and  happens  in  space  uniformly.  In  contrast,  ‘propagating  fragmentation’  designates  a 
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process  wherein  a  medium  gets  fragmented  over  an  extended  duration  of  time  and  in  an 
extended  volume.  The  ‘source  term’  driving  the  fragmentation  is  space-  and  time- 
dependent.  By  maximizing  the  information  entropy  density  at  all  times  and  positions, 
subject  to  constraints  involving  the  source  term,  we  obtain  fragment  size  distributions 
that  are  fully  determined  by  the  source  term  and  the  strain  rate.  We  shall  find  that  weak 
(or  strong)  sources  result  in  a  size  distribution  that  is  sharply  (or  broadly)  peaked  around 
the  most  likely  fragment  size. 

This  is  the  first  instance  in  which  the  duration  of  the  fragmentation  process  enters 
our  consideration.  We  postpone  a  discussion  of  the  temporal  scale  of  fragmentation  to 
a  later  section  (section  4.2)  where  sequential  fragmentation  theories  are  reviewed;  let 
us  note  though  that  mem  can  bear  extension  to  processes  developing  in  time  (Jaynes 
1983,  p  287,  Levine  1985).  The  present  approach  keeps  clear  of  intricate  issues,  such  as 
the  suggestions  that  a  perturbed  system  will  evolve  by  minimal  rate  of  entropy  production 
not  far  from  equilibrium  and  by  maximal  rate  of  entropy  increase  when  far  from 
equilibrium  (Prigogine  1978,  Kaufman  etal  1989). 

We  assume  that  the  fragmented  body  can  be  subdivided  into  regions  of  unit  volume 
that  are  sufficiently  large  for  statistical  considerations  to  be  valid  in  each  region  and  at 
the  same  time  small  enough  that  external  forces  or  stresses  can  be  regarded  as  uniform 
throughout  the  region.  Regions  are  labelled  by  a  convenient  coordinate  r  and  the  time 
by  t.  Fragments  formed  in  (/,  r)  are  characterized  by  linear  size  a  =  a{t,  r)  where  c  is  a 
parameter  specifying  the  fragment.  More  generally  a  might  represent  a  set  of  fragment 
parameters  (size,  shape,  composition,  quality,  etc),  though  we  shall  work  with  a  as 
representing  exclusively  the  linear  size  of  the  fragment.  We  seek  by  mem  the  probability 

p(no  y  ty  r)  with  2  p{na ,  r)  =  1  (31) 

rtfl 

of  finding  n^  fragments  in  unit  volume  around  position  r  and  at  time  t  having  the 
dimension  a.  Under  the  conditions  posited,  the  total  entropy  of  information  is  the  sum 
of  the  entropies  of  the  regions 

S{t)  =  2  S{t,  r)  =  -  2 2 2p(«fl,  t,  r)  logp{n„,  t,  r).  (32) 

r  r  a  tig 

We  maximize  S{t)  for  all  times,  subject  to  constraints  described  in  the  following,  to 
obtain  first  the  probabilities  p  and  then  the  mean  fragment  number  per  unit  volume  for 
fragments  of  type  a, 

n{a,  t,  r)  =  2  ^  0-  (33) 

The  expectation  value  of  any  quantity  Q{a)  is  given  by 

(2>  =  2  Qia)n{a,  t,  r).  (34) 

a 

The  physical  foundations  of  the  method  are  some  theoretical  or  empirical  relation¬ 
ships  between  the  parameter  of  fragmentation  (a)  and  the  parameters  of  the  physical 
process  that  causes  fracture.  The  relations  are  expected  to  have  the  form 

F  (fragment  parameters)  =  (f>  (stress  parameters)  (35) 

of  which  a  particular  case  might  be 

<F(fl,  t))  =  f  « <)  (36) 

where  on  the  left  we  take  an  average  over  a  function  of  fragment  sizes  a  at  time  r  and  on 
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the  right  we  have  a  function  of  the  history  of  the  stress  pulse  crup  to  the  time  t.  Numerous 
forms  of  the  functions  in  (36)  have  been  considered  in  the  past  (Curran  et  al  1987, 
Kachanov  1986). 

We  shall  postulate  an  energy  balance  formulation  (Tuler  and  Butcher  1968,  Kach¬ 
anov  1986)  equating  the  energy  residing  in  the  fragments  with  the  energy  input  invested 
up  to  time  t  by  the  source.  As  already  written  in  (15) 

S  e{a)n{a,  t,  r)  =  E{t,  r)  (37) 

a 

where 

E=WxD  (38) 

=  W^(aI'  [<j((')  -  a,]"  dr-y  (39) 

On  the  left  in  equation  (37)  is  the  sum  over  the  energy  density  functional  for  fragments 
of  size  a,  while  on  the  right  appears  a  semiempirical  relation  for  the  energy  input  in  unit 
volume  required  to  create  the  fragments.  In  (39),  a(0  is  the  stress  at  time  /  and  position 
r,  Oq  is  a  threshold  for  fracture,  (In  a  more  realistic  treatment  one  should  use  several 
thresholds  for  different  kinds  of  stresses:  compressive,  tensile,  shear,)  W  is  the  energy 
at  full  nominal  damage  that  specifies  the  energy  needed  to  produce  the  damage  D  = 
D{t)  defined  by  the  large  parentheses,  A  is  another  empirical  constant  that  depends 
primarily  on  the  properties  of  the  medium:  its  mechanical  strength  and  its  state  (e,g,  the 
population  of  microcracks).  The  power  n  in  the  integral  is  commonly  taken  as  2. 

The  stress  <7(t)  represents  the  external  source  of  fragmentation.  It  is  regarded  in  this 
study  as  an  input  parameter,  being  a  function  of  the  distance  from  the  location  of  the 
energy  deposition  and  of  time.  Theoretical  and  empirical  estimates  of  a{t)  take  into 
account  the  intensity  of  the  source  (e  ,g,  the  powder  factor  of  explosives) ,  the  stress  input 
rate  and  the  elastic  constants  of  the  medium.  The  stress  o  decreases  oscillatorily  with 
both  t  and  r,  due  to  geometry  and  dissipation.  In  truth,  the  energy  loss  due  to  fragment 
formation  itself  reduces  a  and  it  is  necessary  to  use  a  a  that  is  self-consistent  with  both 
the  energy  loss  and  the  degradation  of  the  medium.  For  an  inhomogeneous  medium 
(e,g,  with  substantial  variability  in  crack  densities  inside  it)  A  and  W  will  also  have  a 
distance-  and  time-dependent  behaviour. 

The  energy  functional  e{a)  in  (37)  was  given  in  (10);  we  rewrite  it  as 

e{a)l2jt  =  \pu^a^  +  ha'^a^/K  -  2ya^lR  +  eal2jz  (40) 

where 

=  2jtPa^l5  +  4jtya^.  (41) 

The  separation  of  terms  in  (40)  distinguishes  between  terms  with  dependence  and 
the  two  terms  in  (41)  with  a  different  variation.  Since  by  construction  we  partition  the 
fragmented  material  into  regions  of  unit  volume,  for  which  following  constraint  holds 

{a^}  =  1  (42) 

having  normalized  for  convenience  to  a  cubic  rather  than  spherical  volume,  the  expec¬ 
tation  values  of  the  first  three  terns  in  (40)  are  constant.  Physically  this  means  that  the 
energy  input  for  fragmentation  E  is  fully  expended  on  £„  in  (41),  whereas  the  first  three 
terms  consume  energy  in  a  way  that  is  independent  of  the  state  of  fragmentation. 
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Maximization  of  (32)  subject  to  the  constraints  (37)  and  (42)  leads  to  the  following 
probabilities 

,Ur)  =  exp[~na(^£a  -  va^)]  j 2  exp[-ma(i3e^  -  va^)]  (43) 

where  the  Lagrange  multipliers  and  v  serve  to  satisfy  (37)  and  (42),  The  expected 
number  of  fragments  per  unit  volume  is,  according  to  (33), 

na  =  [exp(i3£„  -  va^)  -  1]"^  f 

The  constraint  equations  give,  for  each  form  of  the  energy  functional  e(a),  a  universally 
valid  solution  as  a  function  of  E  only,  i.e,  independent  of  t  and  r.  The  physical  mechanism 
for  fragmentation  enters  through  the  dependence  of  £  on  r  and  t  (an  example  for  which 
is  the  integral  shown  in  (39))* 

To  describe  the  variation  of  the  Lagrange  parameters  f  and  v  with  £  we  find  it 


convenient  toTtransform  them  to  dimensionless  form, 

j8'  =  Irzya^^  (44) 

v'-a^v  (45) 

where 

a  =  {Sylpe^yl^  (46) 

is  a  length  characteristic  of  fragment  size  (Grady  1982)  shown  in  equation  (26).  With  the 
choice  of  (46)  for  the  length  unit  the  equations  determining  jS'  and  v'  are 

1  =  (s^>  =  2  s^n{s)  (47) 

S 

E'  =  Ejlnya^  =  2  (25^  +  s^)n(s)  (48) 

S 

n(s)  =  {exp[j3'(2s^  +  5^)  -  v's^]  -  1}"’.  (49) 

In  these  units  the  ratio  of  to  a^  is 

e'{s)/s^  =  2s"’  +  s^.  (50) 


As  shown  in  figure  3,  for  v'  tends  to  the  value  -1.0534,  while  varies 

exponentially  as  exp(-l/2l  v'  |£').  the  energy  input  decreases,  increases  and  v' 
changes  sign.  In  the  limit  P'  approaches  infinity  and  v'  approaches  3j3'  from  below. 
However,  the  limiting  value  of  the  energy  £'  where  this  occurs  is  not  zero,  but  3+,  or 
for  the  physical  energy  in  (48) 

E-^En^=6jcya^.  (51) 

This  remarkable  result  of  an  energy  threshold  dictated  by  statistical  requirements, 
and  additional  to  the  crack-mechanical  threshold  implied  by  Cq  in  (39),  will  now  be 
interpreted  on  the  basis  of  equations  (47)-(49).  In  the  simplest  terms,  for  the  mean 
energy  to  be  large,  n  has  to  be  large,  and  this  is  tantamount  to  the  exponent  in  (49)  being 
small.  Therefore  decreases  steadily  with  increase  of  energy  input,  as  seen  in  figure  3. 
However,  the  increase  of  n  has  to  be  kept  in  check,  since  the  expectation  value  of  the 
volume  is  constant  (equation  (42)),  and  this  is  ensured  by  the  other  Lagrange  multiplier 
p'  having  the  tendency  to  counterbalance  the  decrease  of  jS'.  More  illuminating  is  the 
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Figure  3.  Lagrange  multiplier  parameters  in  the  distribution  as  a  function  of  the  energy 
input  £'  responsible  for  the  fragmentation.  The  energy  multiplier  /3'  (left-hand  vertical 
scale)  and  the  volume  multiplier  v'  and  v'  -  3/3'  (right-hand  scale)  are  shown,  all  in 
dimensionless  units  (equations  (44),  (45)  and  (48)).  The  asymptotic  behaviours  of  /3'  and 
v'  are  shown  on  the  far  right. 


behaviour  of  the  distributions  (shown  in  figures  4(fl)-(c))  as  the  energy  E’  varies.  We 
note  that  the  energy  functional  £(.s)  goes  as  and  s^,  while  the  volume  varies  as  Thus 
£(5)  dominates  for  both  large  and  small  sizes ,  while  the  volume  may  do  so  at  intermediate 
ones.  (The  ratio  of  the  two  quantities,  equation  (50),  is  minimum  at  s  =  1.)  Therefore 
the  distribution  is  broad  for  large  energy  inputs  and  sharpens  towards  the  value  s=l 
for  small  E'  (see  figure  4).  At  its  extreme  sharpest,  the  volume  or  energy  distribution 
acts  as  a  delta  function  and  the  ratio  of  the  energy  and  volume  expectations  is  given  by 
3  (equation  (40)  at  j  =  1),  which  represents  the  threshold  for  energy  input  per  volume 
a^.  For  smaller  energy  input  a  maximum-entropy  solution  does  not  exist  and  the  dis¬ 
tribution  is  likely  to  be  an  irregular,  haphazard  one. 

The  general  trend  of  our  results,  including  the  interesting  sharpening  of  the  dis¬ 
tribution  at  small  energy  inputs,  holds  for  energy  functionals  different  from  that  in  (40) 
and  (41),  provided  the  energy/volume  minimum  occurs  at  finite  sizes.  This  permits  the 
generalization  of  our  results  to  include  a  term 

Cs^  (52) 

where  C  is  a  constant  and  2  -  4,5,  which  represents  the  energy  expenditure  in  creating 
internal  damage  (i.e.  cracks  that  do  not  join  up  to  form  fragments)  (Jaeger  etal  1986a, 
Aharony  et  a/ 1986). 

The  distribution  in  (49)  allows,  formally,  fragments  of  size  larger  than  the  unit 
volume,  into  which  the  medium  was  subdivided.  To  lighten  the  paradox,  one  must 
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Figure  4.  Distributions  as  a  function  of  the  reduced  fragment  size  a/a  or  .r;  &  is  shown  in 
equation  (46).  Plotted  are:  (a)  «,  the  expected  number  of  fragments  per  unit  volume 
(equation  (49));  (b)  a^n,  the  expected  volume  (in  units  of  d^);  and  (c)  tafi,  the  expected 
energy  (in  units  of  2jryd^).  In  the  figures:  (i)  is  for  a  strong  source  of  fragmentation,  E'  = 
7.26,  )3'  =  0.001,  v'  =  -1.03;  (ii)  is  for  an  intermediate  strength  source,  E'  =  3.62,  p'  = 
I,  v'  =  2.48;  (iii)  holds  for  a  weak  energy  input  source,  £'  =  3.06,  P'  =  10,  v’  =  29.85. 
The  distributions  sharpen  as  the  energy  input  decreases.  Note  the  changes  in  the  vertical 
scales. 


regard  the  distribution  as  referring  to  an  ensemble  of  fragmenting  bodies,  in  which  the 
occuiTence  of  a  fragment  larger  than  a  unit  adds  to  the  weight  of  distributions  in 
neigbouring  regions. 


4.  Alternative  theories  of  fragment  distributions 

Several  theories  give  n^,  the  distribution  function  of  fragments  of  size  s  (n^  is  the  number 
of  fragments  of  size  s  (^-fragments)).  We  describe  some  of  these.  Summaries  from 
differing  viewpoints  are  due  to  Dehn  (1981),  Grady  and  Kipp  (1985)  and  Campi  (1989). 

4.1.  Geometrical  methods 

4.1.1.  Poisson  statistics.  The  linear,  id  distribution  is  easily  derived  and  will  turn  out 
to  be  of  interest  also  in  higher  dimensions  (in  connection  with  the  so-called  Mott 
distributions).  To  avoid  boundary  effects  we  consider  a  ring  of  perimeter  length  L  which 
is  broken  up  into  P  pieces  by  randomly  placed  cuts  (figure  5).  The  probability  of  obtain¬ 
ing  a  piece  having  a  size  between  s  and  j  +  ds,  anywhere,  is  equal  to  the  product  of 
probabilities  of  placing  one  cut  anywhere  (=1),  of  placing  a  second  cut  in  the  interval 
(s,s  +  dy)  measured  from  the  position  of  the  first  cut  (probability  =  ds/L),  of  not  having 
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Figure  5.  One-dimensional  Poisson  process.  A 
ring  is  broken  up  into  n  (here  9)  segments  by  n 
randomly  placed  cuts. 


Figure  6.  Percolation  construction  of  fragments.  Elementary  segments  of  size  b  are  chosen 
randomly  with  a  probability  p,  resulting  in  compound  fragments  of  sizes  2b,  3b  and  4b. 
An  equivalent  approach  is  to  choose  fragment  voids  with  probability  (1  -  p). 


any  of  the  remaining  P-2  cuts  falling  in  the  interval  (probability  =  [(L  -  s)/LY~'^) 
and  multiplied  by  the  number  of  cuts  (=  P  -  1)  following  the  first.  Altogether 


cbn,(P)  =  1 


ds/L-sV-^ 


(P  -  1). 


(53) 


In  the  limit  of  L/s  and  P  being  large  numbers,  this  reduces  to 

n,  =  (1/5 )  exp(-5/j)  (54) 

where  f=  L/P  is  the  mean  size. 

In  units  where  5  =  1  and  definition  of  a  new  variate  z  through  changes  of  mean,  scale 
and  exponent  according  to 


one  obtains  the  Weibull  distribution 
1  /z-z^yii-' 


exp[- 

where  1/^  is  called  the  shape  parameter. 


Hq\  ) 


-  7  A 


_£o\ 


(55) 


4.1.2,  Percolation  statistics.  In  a  id  discrete  network  consisting  of  L  cells  of  length  b  (as 
shown  in  figure  6)  the  ‘elementary’  segment  sizes  (cells)  are  fixed  {b).  Larger  fragments 
of  size  sb  can  form  by  joinings  contiguous  fragments,  separated  from  neighbouring  ones 
by  the  presence  of  the  edges  of  empty  cells  whose  probability  of  occurrence  on  each  cell 
is  (1  -  p).  The  number  of  s-fragments  is 

n,(p)  =  Lp*-*(l-p)2.  (56) 


For  more  details  on  percolation  in  one,  two,  three  and  higher  dimensions,  see 
Stauffer  (1979, 1985),  Deutscher  etal  (1983)  and  Zallen  (1983). 
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Table  1.  What  is  fixed  and  what  varies  in  Poisson  and  percolation  statistics  (Campi  1988). 
Poisson  Percolation 


Fixed  Number  of  cuts  (?)  Probability  of  segment  (p) 

Variable  Probability  of  segment  Number  of  fragments 


Table  1  sets  out  the  relation  between  percolation  and  Poisson  statistics.  It  can  be 
shown  (Campi  1988)  that  the  mean  or  expectation  number  of  fragments  of  all  sizes, 
namely 

oc 

2  n,(p) 

S=  1 

equals  L(1  -  p)  so  that  in  percolation 

p  =  l-'2n/L  =  l-P/L  (57) 

5 

where  we  have  replaced  In  by  P,  the  number  of  cuts  featuring  in  the  Poisson  statistics. 
Inserting  this  result  in  (56)  we  obtain 

nM^cxp[-{P/L)s]{L/P) 

which  agrees  with  the  result  for  Poisson  statistics.  Size  distribution  of  fragments  (or,  to 
use  the  more  common  term,  of  clusters)  near  the  percolation  limit  p-^  Pc  can  be  shown 
to  be  derivable  from  maximum-entropy  principles  (section  5) .  It  has  also  been  shown  that 
when  the  cells  are  identified  with  crack  positions,  fragmentation  of  a  three-dimensional 
object  occurs,  within  the  percolation  description,  as  a  sudden  transition  with  a  crack 
concentration  p  -  pcii  (Aharony  et  al  1986).  This  is  higher  than  the  usual  critical  con¬ 
centration  for  percolation  pf=  Pci)-  The  reason  is  that  at  pa.  one  gets  a  continuous 
structure  of  cracks  in  the  solid  without  the  solid  breaking  up  into  pieces.  This  (namely, 
fragmentation)  occurs  only  at  a  higher  crack  concentration,  for  p  =  pai  or  higher. 
Computed  fragment  statistics  are  shown  in  figure  7  and  are  compared  with  an  exper¬ 
imental  distribution  (Englman  etal  1984). 

4.1.3.  Tessellation.  Cutting  up  a  plane  with  lines  or  a  solid  with  planes,  or,  alternatively, 
forming  Voronoi  polygons  or  polyhedra,  are  further  ways  of  constructing  fragments. 
The  essential  difference  from  the  physical  standpoint  is  whether  one  allows  lines  or 
planes  to  cut  across  each  other.  (As  a  rule,  true  cracks  do  not  intersect  except  when  their 
propagation  velocity  is  high .)  Extensive  simulations  with  random  lines  in  two  dimensions 
by  Grady  and  Kipp  (1985)  have  indicated  fragment  numbers  that  decrease  exponentially 
either  with  the  linear  dimension  (‘Mott  distribution’)  or  with  the  area  of  fragments.  In 
simulations  by  Sprecher  (Jaeger  et  al  1986a,  c)  the  internal  structure  of  fragments  was 
also  studied.  We  mention  these,  not  only  to  correct  the  imbalance  in  this  review,  which 
deals  almost  exclusively  with  fragment  sizej  whereas  other  parameters  of  fragments  are 
also  susceptible  to  statistical  treatments,  but  also  to  stress  the  importance  of  unfinished 
cracks  within  the  fragment  (‘internal  damage’),  which  can  take  up  a  substantial  part  of 
the  energy  available  for  fracture. 
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than  ^5)  ploited  logantisnically  against 5*-  for  rstenrjtfgts s.  Text  cp^stnm'.td  <^7a 
s£is  are  shou-n  (t^sen  sj-cnbols)  for  the  daal  lattice  {>ercokti(n  jKohanS^  = 0^  03. 
0.4  (so  that  the  crad:  oocupalion  probaMlIlies  are  p  =  O.S.  0.7. 0.6)  ana  data  pi^ts  (fail 
squares)  for  fragnients  larger  than  1  grain,  ohtaioed  froa  the  exploson  of  a  n^tal 
cylindrical  shell  0.63  on  thick  filled  with  47  3DX,/3!  TUTflZ  A!  (Sternberg  1973).  The 
line  designated  ‘Mott  distribution*  represents  equatics  (5S). 


Leaving  aside  the  question  of  physical  justification  for  intersecting  cracks  we  can 
provide  a  heuristic  proof  for  the  Mott  distribution  generalized  to  higher  dimensions 
d  (>1),  namely 

n,  oc  exp[-  (58) 

The  proof  assumes  a  high  density  of  interesting  randomly  thrown  cracks,  such  that  most 
fragments  have  the  minimum  number  of  sides,  i.e.  triangles  in  aplane  and  tetrahedrain 
space.  A  d-sided  object  is  not  yet  a  fragment,  but  must  await  the  arrival  of  another  side 
(the  (d  +  l)th  crack)  to  become  one.  This  will  fall  on  the  d-sided  object  at  a  distance 
from  the  opposite  apex  that  follows  a  Poisson  distribution:  hence  the  distribution  law  in 
(58),  which  is  Poisson-like  in  the  linear  dimension. 


4.2.  Kinetic  distribution 

We  now  describe  theories  in  which  a  sequence  of  fragmentation  stages  is  considered, 
such  that  each  stage  retains  the  memory  of  the  preceding  one  only  by  the  size  of  the 
fragments  that  are  present,  not  through  their  state  of  stress  (McGrady  and  Ziff  1987; 
Cheng  and  Redner  1988,  Derrida  and  Ryvbjerg  1987).  The  results  appear  to  be  inde¬ 
pendent  of  the  physical  causes  of  fracture,  yet  there  is  a  prediction  of  ‘shattering 
transitions’  in  which  the  numbers  of  small  fragments  multiply  to  infinity  as  the  size  of 
the  fragments  decreases. 

The  following  integro-differential  equation  gives  the  number  of  fragments  N{s,  t)  of 
size  s  at  time  t  in  terms  of  the  rate  p{s)  of  5-fragments  breaking  up  and  the  channelling 
ratio  K{s\s)  (this  quantity  gives  the  fraction  of  5'-fragments  (5'  <  5)  created  out  of  the 
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X(s.t)  =  -Hs)V(^.i)-^f  H^ms-,OK(s,sOds-.  (59) 

s 

CbiistT^ts  indcds  mass  conseivation 

r  s'i:(s',s)ds’=s 
-'o 

2sd  2  pt^tulated  mean  aumber  n  offragments  created  at  each  step  (e.g.  2  upon  cutdng) 
OTenby 

r  K(s',s)ds^=/L  (60) 

A  problem  that  can  be  solved  exactly  is 
v(s)^s^ 

(e.g.  A  =  i  for  random  cutting) 

K(s.s’)  =  b(sls’)ls'  (6i) 

(e.g.  6  =  2  for  random  cutting  when  n  =  2  in  (60))- 

The  initial  condition  N(s,  0)  =  d(j  —  L)  for  a  single  piece  has  the  solution 

N(s,  t)  =  e“^[d(j  -  L)  +  2/  -r  f(L  -  s)]. 

The  form  in  (61) 

b(s/s’)  =  (s/sT 

^ves  n  =  (m-T  2)/(m  -}- 1).  For  negative  m  the  number  of  small  jfragments  rises  leading 
to  a  ‘shattering  transition’,  in  which  an  increasing  number  of  tiny  fragments  is  generated 
(Ziff  and  McGrady  1986,  McGrady  and  Ziff  1987). 


4.3.  Physics-based  theories 

To  my  knowledge,  apart  from  the  mem,  which  forms  the  subject  of  this  review,  there 
exist  no  theories  of  fragment  size  distribution  that  are  rooted  in  detailed  physical 
mechanisms.  This  situation  is  probably  due  to  the  difficulty  of  disentangling  the 
stochastic  and  the  deterministic  aspects  of  fracture  (Grady  and  Kipp  1989,  Jaeger  and 
Engiman  1991). 

There  are  indeed  physical  theories  of  fracture,  of  which  Mott’s  is  perhaps  the  most 
notable  (Mott  1947,  IQpp  and  Grady  1985,  Grady  and  Kipp  1989).  These  give,  as  a  rule, 
mean  fragment  sizes,  but  not  distributions. 

Furthermore  there  are  also  models  of  fracture  (mainly  in  brittle  media)  and  of  failure 
(like  yield  in  ductile  solids  (Tate  1967, 1969)).  Some  of  these  are  in  the  form  of  extensive 
computer  codes,  starting  in  the  early  1970s  with  the  one  developed  at  Stanford  (McHugh 
1983,  Seaman  et  al  1984,  Shockey  et  al  1985,  Curran  et  al  1987)  and  simultaneously  at 
Sandia  (Davison  and  Stevens  1972,  Davison  et  al  1977),  and  continuing  with  several 
others  (Rice  1975,  Ravid  and  Bodner  1983,  Ravid  et  al  1987^  Kipp  et  al  1980, 
Taylor  etal  1985,  Kuszmaul  1986,  Brandon  1988).  However,  the  quantitative  measures 
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of  firacture,  failure  or  damage  that  are  grven  by  the  codes,  frequently  by  mapping  them 
onto  the  interior  of  the  disintegrating  body,  are  not  immediately  translatable  into 
fagment  formation  and  into  a  distribution  of  fragment  sizes. 


5.  IVlanmum-entropy  method  and  percolation  theory 


It  is  of  some  theoretical  interest  to  establish  a  link  between  percolation  theory  and  mem, 
though  most  results  are  more  easily  obtained  from  numerical  generation  of  percolation 
clusters  than  by  using  mem  with  the  constraints  imposed  that  are  in  any  case  the  outcome 
of  simulations.  In  an  early  effort  Kikuchi  (1970)  used  a  mean-field  treatment  to  derive 
critical  percolation  probabilities  by  entropy  maximization.  The  agreement  with  correct 
values  was  poor. 

The  events /  in  equation  (1)  for  the  probabilities  p,  are  taken  to  represent  the  finding 
of  rij  clusters  of  size  s  per  lattice  site 

Pj  =P(ns,s). 

We  note  that  in  percolation  a  cluster  arises  as  a  connected  entity  of  bonds  (or  sites)  each 
of  which  is  realized  with  a  fixed  probability  p.  The  mean  number  of  s-clusters  is 


<«,>=  E  p{ns,s)nj. 

n,=0 

The  constraints  are,  as  before, 

Ep(«5,5)  =  1 

for  eachs  and  the  following  relations  are  given  by  Stauffer  (1981)  to  define  critical  indices 
or,  7  and  6: 


(  S  <«,)) 

>  Sing 

(62) 

'  j  '  sing 

(63) 

s 

(64) 

E  5(l-e"'‘")nj(p-j-Pc)oc/ii/«, 

(65) 

In  these  relations  e  =  [p  -  pd,  where  pc  is  the  critical  occupation  probability  for  per¬ 
colation,  and  the  subscript  ‘sing’  represents  the  leading  singular  part  of  the  sum  which 
remains  finite  at  p  =  p^. 

The  MEM  solution  for  the  cluster  size  distribution  is 
p(/2„  j)  =  L-^  exp(-/7,x,  -  A‘) 

where 


Xs  =ft-»o  A  +  Bs+  Cs^  +  7)5(1  - 

(fij)  =  -  1)“‘ 


(66) 


and  A'  is  a  normalizing  constant  independent  of  the  size  L.  The  coefficients  -A, . . 


D 
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are  fitted  to  satisfy  the  constraints  (62)-(65)-  The  fitting  is  not  easy  to  carry  out  since 
several  limiting  processes  are  involved,  namely  £->0,  ^->0,  L 

(size  of  lattice)— >  however,  it  is  clear  that  for  5  in  some  region  {s  large,  yet  smaller 
than  L^)  the  terms 

Cs-  -r  Ds(l  -  e-*") 

will  dominate  the  expression  (66)for  Rewriting  this  for  the  given  range  of  s  as 

jc,  Cs^  +  ^  Oih)  (67) 

and  noting  that,  for  x,  1, 

<«,)  =  L-^x;^  (68) 

we  conclude  that  the  mean  number  of  clusters  of  size  s  will  show  an  effective  power-law 
dependence  of  the  form 

{ns)^s-^  (69) 

where  2  <  t  <  3. 

For  another,  heuristic  derivation  of  this  result  see  Zallen  (1983).  Numerical  esti¬ 
mation  of  cluster  numbers  near  p  =  Pc  gives  in  three  dimensions  t  =  2.2,  whereas  more 
generally  (with  p  p^)  the  relationship  in  (69)  becomes  modified  by  a  smooth  bell-like 
function  of  s  (Stauffer  1979).  However,  outsize  cluster  numbers  decrease  with  size 
exponentially. 

The  relation  betw'een  fragment  sizes  and  percolation  clusters  appears  to  have  been 
determined  in  the  first  instance  by  Englman  et  al  (1984).  A  fragment  is  formed  when  a 
piece  of  material  is  fully  surrounded  by  linear  cracks  (in  two  dimensions)  or  by  crack 
plaquettes  (in  three  dimensions).  The  size  of  fragments  is  given  by  the  cluster  size  in  the 
dual  lattice,  defined  by  placing  a  bond  across  a  face  where  a  plaquette  is  missing  (or 
across  a  line  where  a  crack  is  missing,  in  two  dimensions)  and  removing  a  bond  where  a 
crack  plaquette  is  present.  The  corresponding  occupation  probabilities  in  the  two  lattices 
are  related  by 

Pdual  1  Paack* 

Fragment  distributions  for p„ack  >  Pc  (=  the  critical  crack  concentration)  were  obtained 
by  a  two-dimensional  simulation  and  followed  a  modified  Mott  distribution  law,  with 
positive  or  negative  deviations  at  very  large  sizes  (Englman  et  al  1984)  (figure  7). 

Let  it  be  remarked  that  mechanical  failure  of  the  solid  occurs  at  p^ack  =  Pcii>  so  that 
the  condition  Pc^ack  ^  Pen  will  be  achieved  in  practice  through  continued  stress  loading 
of  the  solid,  while  it  is  being  held  together  artificially:  either  by  confinement  or  inertially. 
In  an  impact  loading  (by  projectiles,  fragments  or  laser  light)  the  confinement  would 
not  be  expected  to  occur  and  the  fragmentation  would  take  place  as  soon  as  p^ack  =  Pcii 
obtains  (Blackman  and  Goldsmith  1978,  MacAulay  1987).  Then  one  would  expect  to 
see  the  relation  (69)  in  effect. 


6.  Experimental  distributions  in  relation  to  mem  predictions 
6.1.  Exploding  cylindrical  shells 

Cylinders  made  of  Armco  iron  or  of  heat-treated  steel  were  filled  with  explosives  and 
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Figure  8.  Cumulative  fragment  distributions  from  two  tests  by  Mock  and  Holt  (1983)  on 
metal  cylinders.  Total  distributions  are  shown  together  with  partial  distributions  for  types 
1-4,  classified  according  to  fragment  shapes  and  origin.  The  broken  curves  are  theoretical 
fits  in  which  the  contribution  of  the  strain  energy  term  (equation  (10))  changes  with  the 
type  number  as  follows:  1  2  3.  The  full  curve  for  type  4  guides  the  eye. 


the  resulting  fragment  size  distributions  analysed  by  Mock  and  Holt  (1983)  (figure  8). 
The  cylinder  walls  were  about  2  cm  thick  and  the  mean  diameter  about  20  cm,  which 
makes  it  questionable  whether  to  apply  two-dimensional  (shell-like)  or  three-dimen¬ 
sional  modelling.  However,  the  authors  separated  the  fragments  in  accordance  with  a 
characterization  scheme  based  on  the  nature  (brittle  or  shear)  of  the  fracture  surfaces. 
Different  fragment  types  (labelled  type  1-4  in  figure  8)  had  different  mean  mass  per 
fragment.  The  type  with  the  largest  mean  mass  (type  1  in  figure  8)  followed  a  roughly 
exponential  cumulative  distribution,  except  for  small  sizes,  as  a  function  of  the  fragment 
mass  M.  However,  the  types  with  smaller  mean  masses  (types  2-4)  showed  a  decrease 
with  mass  that  was  faster  than  that.  The  explanation  is  straightforward  in  terms  of  the 
formulae  (10),  (18)  and  (19).  The  small  mean  mass  is  indicative,  by  (26),  of  larger  strain 
rate  e.  In  formula  (10)  for  the  exponent  in  the  distribution  (/lo),  equation  (19) ,  the  strain- 
rate  term  goes  as  a^,  or  the  5/3  power  of  mass,  to  be  compared  with  the  inertial  term 
going  as  the  mass,  which  is  expected  to  be  dominant  for  the  smaller  strain-rate,  type  1, 
fragments. 

A  large  number  of  fragment  distributions  arising  from  explosively  shattered  cyl¬ 
indrical  shells  were  obtained  by  Sternberg  (1973)  with  variations  in  the  shell  composition 
and  treatment  and  in  the  explosive.  A  simple  exponential  (Mott-type)  distribution  failed 
to  account  for  the  data  over  the  whole  mass  range  and  consequently  a  distribution  was 
synthesized  from  three  separate  exponentials,  appropriate  to  small,  medium  and  large 
sizes.  The  mem  distribution  in  two  dimensions  is  given  as  a  function  of  linear  dimension 
a  by  equation  (19),  where  Xg  takes  the  form  analogous  to  equations  (10)  and  (18) 

Xa^Aa-¥  Ba^  4-  Ca^.  (71) 

With  suitable  choice  of  the  parameters  A,  B  and  C  the  mem  fits  quite  well  the  results  of 
Sternberg  (1973)  for  two  combinations  of  steel  and  explosive  (figures  9(fl)  and  (b)). 

Exploding  mild-steel  tubes  with  wall  thickness/diameter  ratios  ranging  from  0.05  to 
0. 17  were  investigated  by  Stronge  et  al  (1989) ,  regarding  the  dependence  of  the  fragment 
size  distribution  on  the  fracturing  mechanism.  This  was  varied  through  changing  the 
explosive  (charge)  mass/metal  mass  (C/A/)  ratio  over  a  10-fold  range.  A  two-dimen- 
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Figure  9.  Distributions  of  fragments  from  an  exploding  metal  shell  versus  the  square  root 
of  the  fragment  mass  (Sternberg  1973).  The  broken  curves  show  two-dimensional  mem 
distributions,  equation  (71),  in  which  represents  the  fragment  mass  and  the  parameters 
take  the  following  values  (in  mass  units  of  milligrains):  (a)  A  =  0.22,  B  =  0.05,  C  =  0.05; 
(b)  A  =  0.22,  B  =  0.05,  C=  0.0625. 
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Figure  10.  Mean  mass  m  of  mild-steel  fragments 
versus  the  charge/metal  mass  C/M  ratio  (Stronge 
et  al  1989).  The  broken  curve  siiows  the  theor¬ 
etical  result  of  equation  (73). 


sional  Mott  distribution  was  assumed  without  obvious  evidence.  The  mean  mass 
decreased  with  C/M  as  seen  in  figure  10. 

In  the  theory  the  most  likely  mass  m{^a^)  varies  with  the  strain  rate  as 

m  oc  £-2.  (72) 

We  can  relate  the  explosive/metal  mass  ratio  to  the  strain  rate  by  using  the  Gurney 
relation  (Jones  etal  1980): 

Uo  =  [2G/(0.5  +  M/C)]‘/2  (73) 

where  Uq  is  the  velocity  with  which  the  outer  part  of  the  cylinder  moves  and  G  is  the 
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energy  release  by  the  explosive  per  unit  mass.  We  associate  the  strain  rate  with  VQ/innev 
cylinder  radius  and  obtain  from  (72)  and  the  expression  for  Vq 

m  oc  0.5  + {MIC).  ■  (74) 

The  broken  curve  in  figure  10  indicates  quite  reasonable  agreement  with  the  exper¬ 
imental  data  points. 

On  the  basis  of  a  large  number  of  fragmentation  experiments  with  exploding  cylin¬ 
ders,  Lin  et  al  (1989)  proposed  a  numerical  fit  for  the  number  of  fragments  as  a 
function  of  the  dimension,  hardness  and  composition  of  the  cylinder  and  of  C/M,  the 
charge/metal  mass  ratio: 

=  lilt  +  107CP  -  0.39H  +  (C/M)(-3250  +  10462CP  -  2.52H)  +  34.5 

where  t  is  the  thickness  of  the  cylinder  (cm) ,  H  its  (Vickers)  hardness  and  CP  the  carbon 
content  percentage,  (It  is  appropriate  to  remark  here  that  Mott’s  involvement  in  the 
mechanism  of  fracture  arose  during  the  Second  World  War  from  his  investigation  into 
the  difference  between  iV,ot  in  Allied  and  German  shells.  Eventually  it  was  all  traced  to 
CP,  which  increased  the  metal  brittleness)  (Mott  1984,  private  communication),) 


6.2.  Fragments  in  mining  operations 

It  might  be  thought  that  a  distribution  with  three  free  parameters,  such  as  given  by 
equations  (10)  and  (19),  is  bound  to  agree  with  practically  any  experimental  curve, 
especially  when  the  cumulative  distribution  is  plotted  (on  a  logarithmic  scale).  The 
following  example  shows  that  this  is  not  the  case  and  that  when  sufficiently  accurate  data 
are  plotted  for  a  wide  enough  fragment  size  range  the  mem  distribution  postulated  in 
equation  (10)  can  be  discredited  and,  instead,  another  mem  distribution  with  different 
constraints  appears  to  operate.  Unfortunately,  the  latter  constraints  do  not  have  any 
simple  physical  meaning. 

Boreholes  of  3.2  cm  diameter  were  filled  with  Donarite  and  exploded  in  an  under¬ 
ground  site  of  the  Rotem  (Negev,  Israel)  oil-shale  fields  of  the  PaMA  company  (Jaeger 
et  al  1986b).  Post-explosion  fragments  were  kept  in  place  by  chicken  wire,  collected 
from  several  locations  and  layers  around  the  original  boreholes,  the  pieces  individually 
weighed  and  their  distributions  determined  for  each  location.  The  derived  distributions 
differed  significantly,  according  to  the  distance  of  the  locations  from  the  centre  and 
somewhat  on  their  orientations  (the  terrain  was  anisotropic  due  to  major  faults). 
However,  (quite  remarkably)  all  distributions  could  be  brought  into  coincidence  by 
scaling  the  sizes  at  each  location  with  a  scale  factor  X{R)  that  increased  monotonically 
with  the  distance  R  of  the  location  from  the  borehole  as  (figure  j j, 

We  have  not  found  it  possible  to  fit  the  collapsed  distribution  with  any  set  of 
Lagrangian  parameters  appearing  in  equation  (19).  A  least-squares  fit  of  the  parameters 
was  not  done,  but  the  intuitively  ‘best’  fits  shown  in  the  figure  as  broken  curves  are 
evidently  inadequate  to  reproduce  the  results  over  the  four  decades  of  linear  size 
variation.  A  MEM  distribution  of  the  form 

n,  =  {exp((a/a,)  +  b  log(fl/flo)]  -  1}"‘  (75) 

was  found  to  fit  the  observed  cumulative  weight  distribution 

j  pd^ng  da 

reasonably  well,  with  the  choice  and  interpretation  of  the  parameters  as  follows: 


123 


Fragments  of  matter  from  a  maximum-entropy  viewpoint 


Figure  II.  Cumulative  weights  of  mined  oil  shale  against  normalized  linear  fragment  size 
(Jaeger  el  al  1986b).  Results  from  five  regions,  numbered  according  to  increasing  distance 
from  the  borehole,  He  on  the  same  curve  when  a  length  scale  parameter  K{R)  is  adjusted. 
The  broken  curves  are  based  on  equations  (18)  and  (19),  and  are  plainly  inadequate.  The 
full  curve  is  obtained  using  equation  (75)  with  parameters  given  in  the  text. 

flj  (upper  cut-off  of  size)  =  2K{R) 

Uq  (lower  cut-off  of  size)  =  8  x  10"^A(/?) 
b  (effective  inverse  power  of  distribution)  =  3.43. 

The  two  functions  in  (a  and  log  a)  indicate  the  existence  of  physical  constraints 
(appendix  1)  on  the  distribution.  The  constraint  on  a  can  be  accorded  some  physical 
meaning  in  terms  of  dense  cracks  (along  the  lines  of  the  discussion  in  section  4.1.3); 
however,  the  constraint  related  to  log  a  is  obscure. 

In  a  similar  type  of  experiment,  though  on  a  much  smaller  scale,  Curran  et  al  (1977) 
collected  the  ejecta  from  the  crater  formed  by  high  explosives  in  contact  with  the  surface 
of  fine-grained  quartzite  rock.  The  number  of  fragments  per  unit  fragment  radius  a,  n^, 
follows  (for  about  five  decades  of  variation  in  nj  remarkably  well  the  power-law 
relationship 

n^  a  *  (76) 

indicating  the  importance  of  the  surface  energy  term  (proportional  to  a^)  in  equation 
(1),  (Figure  4.10  in  Curran  et  al  (1987)  shows  the  cumulative  number  of  fragments. 
The  theoretical  fit  in  the  figure  suffers  from  identification  of  fragment  and  crack  size 
distributions  in  their  equation  (4.6).)  The  surface  term  is  indeed  expected  to  dominate 
for  small  fragments  (a  <  a)  and  for  such  fragments  we  anticipate  a  power-law  behaviour, 
like  that  in  (76). 

6.3.  Distribution  of  space  debris 

Man-made  debris  in  space  exceeds  the  number  of  natural  meteoroids  and  represents  a 
danger  to  spacecraft  (Kessler  1985,  Rajendran  and  Elfer  1989).  The  fragments  arise 
from  defunct  space  vehicles  that  have  broken  up,  exploded  or  become  involved  in 
collisions.  Added  to  this  are  solid  fuel  particles.  Three  distributions  are  given  as  a 
function  of  the  mass  m  (Johnson  and  McKnight  1987)  or  of  the  size  s. 
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Figure  12.  Distribution  of  masses  in  space  due 
to  defunct,  broken-up  satellites.  Broken  curve: 
observation  (described  in  the  text).  Full  curve: 
fit  by  a  power  law  based  on  percolation  theory. 


Small  fragments  due  to  collisions  have  a  power-law  (fractal-like)  distribution. 

n,  =  (77) 

Larger  fragments,  but  not  exceeding  1960  g,  follow  the  distribution 

=  C2fn~^^^  exp(-0.0575w'/^)  (78) 

with  m  in  grams.  For  even  larger  fragments,  caused  mainly  by  explosions,  one  finds 

n,  -  exp(-0. 0205m '/^).  (79) 

The  coefficients  Cj,  C2  and  C3  depend  on  the  total  masses  included  in  the  distribution 

and  are  not  relevant  to  the  expected  distribution  from  a  single  event.  In  figure  12  we 
have  combined  the  three  distributions  in  one  smooth  curve  and  compared  the  resultant 
with  the  percolation  result  n, « 

The  experimental  distribution  drops  slower  than  this  prediction  for  small  sizes  and 
faster  for  larger  fragments.  In  some  sense,  the  prediction  represents  a  simple  averaging 
over  the  full  range. 

6.4.  Charged  atomic  clusters 

Clusters  consisting  of  10-10^  atoms  or  molecules  are  of  interest  theoretically  (as  tran¬ 
sition  species  between  molecules  and  solids)  and  experimentally  (e.g.  as  potential 
catalysts).  A  review  of  early  results  was  given  by  Phillips  (1986).  Clusters  can  be 
formed  from  the  gaseous  phase  by  accretion  or  condensation  after  passing  through  a  jet 
(Levinger  et  al  1988 ,  Rayane  et  al  1989)  or  from  an  impacted  solid  by  sputtering  ( Weiland 
et  al  1989).  The  distribution  of  cluster  masses  can  be  obtained  by  time-of-flight  mass 
spectrometry  and  depends  on  the  conditions  under  which  clusters  form  (e.g.  the  ambient 
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gas  pressure).  The  distributions  are  frequently  peaked  (as  in  rare-gas  clusters  Ar^)  at 
magic  numbers  associated  with  favourable  conditions  for  packing  or  with  bunching  of 
electronic  (binding)  energies  (Knight  etal  1984,  Stephens  and  King  1981).  The  broad 
envelopes  of  mass  distribution,  \  hich  are  our  concern  here,  generally  possess  a 
maximum,  though  exceptions  from  this  situation  do  exist.  Thus,  counts  of  positively 
charged  uranyl  acetate  clusters  decrease  monotonically  for  cluster  sizes  of  up  to  50  units , 
while  the  negatively  charged  clusters  show  a  maximum.  Charged  clusters  are  easier  to 
observe  in  experiments  than  neutral  ones. 

Writing  the  mem  distribution  as  a  function  of  the  mass  m  in  the  form  appearing  in 
equation  (19),  namely 

„^  =  (e.(m)_i)-l  (gO) 

where 

x{m)  =  X  drtti'  (81) 

r  " 

(cf  equation  (19)),  we  note  that,  for  n^  to  have  a  maximum  at  a  non-zero  mass,  at  least 
some  of  the  exponents  in  the  above  sum  must  be  negative.  We  recall  that  the  energy 
equipartition  constraint  of  Griffith  (1943)  described  in  section  3.1  led  to  a  negative 
exponent  r=  -1/3. 

This  arose  from  supposing  that  the  available  energy  was  first  distributed  equally 
between  all  atoms  or  molecules,  so  that  each  atom  absorbs  a  part  of  the  order  m~*  and 
then  expends  this  on  the  surface  energy  (which  is  proportional  to  In  addition,  one 
can  suppose  that  the  absorbed  energy  is  also  expended  on  creating  a  transverse  dynamic 
tension  (analogous  to  the  strain  energy  term  in  (10))  proportional  to  Then  in  the 
sum  (81)  one  also  has  a  term  with  exponent  r  =  2/3. 

A  fitting  of  expres.sions  (80)  and  (81)  to  the  mass  distribution  data  of  Levinger  et  al 
(1988)  for  Ar^  was  made,  with  the  choice  of  the  coefficients  in  (81),  given  (in  Ar  mass 
units)  by 

fl-l/3  =  8.5  02/3  ==  15 

(figure  13).  Fits  of  similar  goodness  were  achieved  to  other  data  obtained  by  the  above 
authors  for  Ar^  and  by  Rayane  et  al  (1989)  for  positively  charged  indium  and  lead 
clusters,  by  choice  of  slightly  different  coefficients.  The  formation  of  small  ions  (m  <  5) 
is  probably  due  to  factors  not  included  in  the  energy  constraint  and  for  these  the  misfit 
is  excusable .  It  is  notable  that  the  skewness  of  the  experimental  distribution  is  reproduced 
by  the  theory.  To  end  this  section,  we  wish  to  stress  this  adaptation  of  the  mem  to 
fragment  formation  by  growth,  rather  than  by  disintegration,  as  in  the  rest  of  this  article. 


d.5.  Droplets 

The  distribution  of  droplets  in  a  spray  is  important  for  the  operation  of  fuel  injection 
engines,  for  painting  and  elsewhere.  The  liquid  column  or  sheet  is  unstable  due  to 
fluctuations  and  breaks  into  droplets  that  are  subject  to  several  forces  (inertial,  nozzle 
pressure ,  aerodynamic  drag  by  the  ambient  gas) .  Studies  have  been  made  of  droplet  sizes 
and  velocities  as  a  function  of  pressure  and  distance  along  the  flow  and  perpendicular  to 
It.  Experimental  investigations  use  phase/Doppler  techniques  (Bachalo  et  al  1988, 
Sellens  1989).  For  a  description  of  the  subject,  including  an  introduction  to  size  dis¬ 
tribution  results,  the  book  of  Lefebvre  (1989)  is  of  use.  mem  were  applied  to  obtain  the 
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Figure  13.  Counts  of  charged  argon  clusters  against  argon  mass  number  n.  Full  circles; 
experimental  data  obtained  by  time  of  flight  (Levinger  et  al  1988).  Open  circles:  theoretical 
results  (equations  (80)  and  (81)  with  parameters  shown  in  text). 


Drop  size  (uin) 

Figure  14.  A  typical  surface  of  droplet  distributions  as  a  function  of  linear  size  and 
downstream  velocity.  (After  Sellens  (1989).) 


joint,  velocity  size  distribution  function,  employing  constraints  on  the  total  mass,  the 
surface  energy,  kinetic  energy  and  total  momentum  (Sellens  and  Brzustowski  1986, 
Sellens  1989).  It  is,  however,  stated  that  ‘it  is  simpler  to  measure  the  resulting  spray, 
rather  than  the  input  conditions’  (Sellens  1989).  An  example  of  the  experimentally 
determined  joint  distribution  is  shown  in  figure  14  as  a  function  of  the  droplet  size  6 
(typically  10^  lim)  and  velocity  y.  The  theoretical  distribution  has  the  form  shown  in 
equation  (A1.3),  where  0  are  functions  of  d  and  v  entering  the  constraint  expressions. 
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The  marginal  size  distribution  reproduces  the  observations  satisfactorily,  but  in  the 
v-d  correlations  there  are  discrepancies.  The  use  of  correlations  in  mem  has  been 
described  (Englman  et  al  1988b).  In  the  experiments  the  mean  velocity  decreases  with 
droplet  size,  due  to  the  aerodynamic  drag  being  more  effective  for  small  droplets,  while 
the  prediction  shows  constancy  (since  the  drag  is  not  formulated  as  a  constraint). 
This  situation  differs  from  that  typical  for  fragmentation  of  solids,  where  smaller-sized 
fragments  fly  off  faster  than  larger  ones  (Grady  et  al  1985).  The  latter  result  can  in  fact 
be  regarded  as  a  direct  consequence  of  the  momentum  and  kinetic  energy  constraints, 
which  introduce  the  functions 

0  =  vd^  and 

in  the  exponential  in  ( Al .  3) ,  so  that  (the  Lagrange  multiplier  A''  being  presumed  positive) 
smaller  sizes  lead  to  higher  speeds. 

In  a  related  subject,  although  many  orders  of  sizes  away,  the  condensation  of 
interstellar  clouds  by  fluctuations  was  considered  (Kiguchi  et  al  1988,  Narita  etal  1988, 
Stabler  1983).  Blobs,  stars  of  10“^-10"*  solar  masses,  materialize  by  collapse  and  frag¬ 
mentation,  and  so  do  rings  of  up  to  10^  solar  masses  from  rotating  clouds. 

Returning  to  Earth,  we  note  the  exponential  distributions  with  a  low  size  cut-off  in 
one-dimensional  composite  fibres  (Wagner  and  Eitan  1990). 

6.6.  Nuclear  disintegration 

A  percolational  critical  point  has  been  shown  to  operate  in  nuclear  disintegration 
(Campi  1988,  Gross  1990).  High-velocity  gold  ions  impinging  on  emulsion  produced  a 
distribution  of  fragments  with  charges  between  1  and  79  (=  the  number  of  charges  on 
Au),  observed  by  Waddington  and  Freier  (1985).  Campi  (1988)  has  found  correlations 
between  the  second,  first  and  zeroth  moments  of  the  size  distribution  that  strongly 
resembled  the  correlations  in  a  percolative  system  near  p  =  p^.  (Use  was  made  of  the 
relation  (57)  originally  noted  by  him.) 

The  largest  fragment  in  disintegration  can  be  associated  with  the  cluster  that  tends 
to  percolate  (i.e.  extend  throughout  the  media)  at  critical  concentration  p^  and  whose 
mass  density  is  known  as  the  fractal  dimension  (Zallen  1983).  In  a  finite-sized  nucleus 
the  largest  fragment  size  is  related  to  the  total  mass  of  the  nucleus  by  finite  size  scaling 
(Stauffer  1979).  By  proper  identification  of  the  largest  fragments  and  of  p^  in  the 
experimental  data  of  several  atoms,  Campi  (1989)  has  found  a  fractal  dimension  that 
agrees  with  that  for  the  percolating  cluster.  In  the  same  paper  Campi  proposed  a  valuable 
hypothesis,  elaborated  by  Gross  (1990),  for  the  ergodicity  of  the  percolation  model, 
which  is  also  pertinent  to  the  ergodic  behaviour  underlying  the  use  of  mem.  According 
to  this,  the  sudden  fracturing  event  leaves  the  incipient  fragments  in  an  excited  state, 
which  probes  all  distributions  in  an  ergodic  fashion.  At  a  later  stage,  a  relaxation  occurs 
in  which  mass  distributions  do  not  change  any  more. 


7.  Conclusions 

Fragment  distributions  based  on  maximum  entropy  not  only  have  theoretical  soundness 
(alas,  not  universally  acknowledged) ,  but  they  also  possess  the  merits  of  containing  both 
probabilistic  and  physical  components  and  of  accommodating  naturally  exponential 
and  inverse  power  behaviours  (as  a  function  of  particle  size).  Various  experimental 
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distributions  have  been  well  reproduced  by  mem,  though  other  ad  hoc  laws  can  do  this 
too,  except  in  the  extremities  of  very  small  and  large  fragments.  Some  characteristics  of 
empirical  distributions  (like  fractal  nature  or  possession  of  a  hump)  have  their  signature 
in  MEM,  leading  (with  some  measure  of  certainty)  to  physically  significant  conclusions 
as  to  the  breaking  mechanism. 

Concerning  useful  extensions  of  the  subject,  one  envisages  future  uses  of  maximum 
entropy  for  non-equilibrium  phenomena,  like  rates  of  catalytic  reactions  involving 
particulates  or  fluid  flow  in  fragmented  media. 
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Appendix  1.  On  maximum  entropy 
ALL  Basics 

In  the  basic  formulation  of  the  maximum-entropy  method  (mem)  one  maximizes  the 
information  entropy  S  with  respect  to  the  probabilities  p,  of  events  /,  subject  to  the  R 
observational  inputs  or  constraints  written  as: 

<00  =  2p/0'’(/)  =  f  p{j)VU)  d/  r  =  1, 2, . . . ,  R.  (Al.l) 
/  •' 

Summation  or  integration  (for  continuous  variates)  is  over  the  range  of  variation. 

Formally,  and  in  practice,  one  maximizes  the  following  function  of  probabilities, 
named  Lagrangian  (Rivier  etal  1989): 

Uip})  =  -^Pi  logR/  -  ^^i^Pj4>\j)  -  <0‘>) 

-  ...  -  X^(^Pj<p^{j)  -  (0'^))  (A1.2) 

obtaining  for  the  mem  probabilities 

Pj  =  exp[-A‘0H/)  -  ...  -  A^0^(/)]  (A1.3) 

in  which  the  Lagrange  multipliers  are  fitted  to  satisfy  the  constraints  (Al.l).  As  already 
noted  in  the  introduction,  the  mem  solution  is  characterized  by  a  smooth  behaviour  as  a 
function  of  j. 
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Examples  of  constraints  are  first,  the  normalization  condition  on  the  probabilities, 
l  =  2p/=  [/?(/■)  dy  (A1.4) 

and  observed  averages  like 

{})  =  2  Pjj  and  (j^)  =  S  Pjj'^ 

i  ! 

etc.  Acceptable  constraints  are  also  observed  frequencies  of  events,  in  the  form 

(Al,5) 

i 

Because  of  the  presence  of  the  singular  function  djj^  or  6{j  -  jo)  in  the  exponent  of  pj  in 
(A1.3),  the  MEM  solution  is  in  general  no  longer  smooth,  but  is  probably  still  more 
conservative  than  alternative  non-MEM  solutions. 

In  fact,  the  range  of  the  variate  j  can  also  be  written  as  a  constraint,  in  the  sense  of 
zero  frequency  of  occurrence  outside  the  range. 


AL2.  Priors  n(j) 

These  are  introduced  as  generalizations  of  the  summations  over  the  events  in  (A1.2), 
with  the  purpose  of  weighting  each  event  by  a  factor  n(y),  so  that  the  entropy  becomes 

5n  =  -  S  n(y)/7;  \ogpi  (A1.6) 

/ 

and  a  typical  constraint 

{<t>)n  =  2  n(j}pi<l>{i).  (A1.7) 

The  priors  are  weights  ascribed  a  priori  (that  is,  without  regard  to  observations)  to 
each  event  and  correspond  to  ‘measure’  in  set  theory.  Skilling  (1989)  writes  about  them 
thus:  ‘n  is  the  Lebesgue  measure  associated  with  [the  variable],  which  must  be  given 
before  an  integral  can  be  defined.’  A  possible  probabilistic  interpretation  identifies  the 
priors  with  subjective  beliefs  or  reliance  attached  to  the  events.  Physically,  priors  may 
be  related  to  instrumental  effects.  An  example  of  this  is  the  photofragmentation  studies 
of  Ar^  cluster  ions  (section  6.4).  In  making  a  statistic  of  the  number  of  detached  ions  as 
a  result  of  photoabsorption,  the  data  must  be  weighted  by  the  dependence  of  the 
photoabsorption  cross  section  on  the  size  (n)  of  the  cluster  ion  (Levinger  et  al  1988). 

If  the  priors  are  known  or  postulated  there  is  no  difficulty  since  one  can  proceed  as 
follows.  Let  the  weighted  probabilities  of  events  be  denoted  by 

Pj  =  nif)pf.  (A1.8) 

Then  we  rewrite  (A1.6)  and  (AL7)  as 

5n  =  -2i’,logP;/nO-)  and  Wn  =  2  i’/d'O).  (A1.9) 

J 

We  maximize  the  Lagrangian  with  respect  to  P,.  Equation  (A1.8)  is  the  solution, 
with p,  given  in  (Al  .3).  The  functional  form  of  the  exponential  part,  i.e.  the  dependence 
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oa  o^,  g\ . . is  sot  sgected  fay  ifee  istrodiar^oa  of  ifas  piicr,  fact  tlie  pie- 
expoasntialpaitis,  aad3oaielfaeE22Esn2alva2!i£SofliisLagiassa3iasl!^!:e4S,  A*  to 
7A  This  follows  ^ce  these  aie  ^'ea  the  ia^lidt  eqosdoas 

{&)  =  S  e3q)I-A  V‘0)  -----  VOll  (-A1-10) 
/ 

which  contsm  II,  In  the  use  of  B.JEi5  as  a  pare  predictive  tool  for  natural  phenomena,  the 
dependence  of  the  soludoa  oa  the  priors  raises  the  question:  ‘Wiudi  is  the  best  prior?’ 
(JeSreys  1961,  Jaynes  1983,  Enghnan  et  al  1988a).  It  seems  that  no  saJisfeatory  answer 
exists  at  this  stage,  except  perhaps  to  say  that  the  solution  is  unique  as  ^  as  its 
asymptotic  exponential  dependence  is  coacemed  but  not  regardingthe  piefactor. 

The  following,  related  problem  appears  to  have  been  solved  piivier  et  ed  1989). 
Suppose  that  there  ejdst  a  number  K  of  pos^Ie  priors  (e.g.  K  alternative  belies  or 
hypotheses),  given  by 


n^O  k=h....K. 

(A1.U) 

What  linear  combination 

K 

n(j)  =  2  p‘n»a) 

k=l 

(A1.12) 

gves  the  best  prior?  ‘Best*  is  to  be  considered  in  terms  of  incoming  information,  since 
before  that  information  we  suppose  (democratically)  that  all  priors  have  the  same 
weight.  Theproposedsolutionreliesonmaximidngthefollowinggeneralized  functional 
of  the  event  probabilities  and  of  the  weights 

^({P/h  {?'}}  =  -  2p/  log  Pi/no)  -  2  ;.'(2  P/^'U)  -  m)  - 

A 

r/  2  p*  logp^ 

1 

(A1.13) 

where  II  is  defined  in  (A1.12)  and  v  quantifies  the  credit  accorded  to  the  diversity  of 
views  relative  to  factual  evidence.  (Rivier  et  al  (1989)  assumed  =  1,  which  is  as  good 
a  starling  point  as  any.)  The  MEM  weights  are 

=  exp{Xklv)/'^e.xp{Xjv) 

'  k 

(A1.14) 

where 

Xi,  =  2  TL^(J)  exp(-  2 

(A1.15) 

is  the  overlap  between  the  prior  II*'  and  the  exponential  part  of  the  mem  probability.  The 
dependence  of  the  weights p’‘  on  the  incoming  evidence  is  through  the  expectation  values 
{(j)’)  and  the  resulting  values  of 


It  is  interesting  to  compare  this  result  with  the  Bayesian  rule  (Jaynes  1985,  Jeffreys 
1961,  Skilling  1989) 

P(k\{(P))  =  P{k)P{mk)/Pm  (A1.16) 

whose  meaning  in  the  present  context  is  the  following.  The  probability  P{k\{<p))  for  the 
hypothesis  (or  prior)  k,  given  that  the  mean  has  the  observed  value  (0),  is  equal  to 
the  probability  P{k)  of  the  same  hypothesis,  irrespective  of  the  observation  on  {(p). 
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Eialtipiisd  bj’  ifee  protebfSly  P{ipy.  kj  of  cblsssing.  fp)  ca  tbs.  ba^  of  tbs  k  bypolbs^ 
sad  dhidsd  by  tba  prcbabSi^’  of  oblsiaiag  lire  obseivaiioa  (p}  on  tbe  nieaa.  Tbs 

Bayssssa  ialsaprstsifoa  of  ibs  rssaii  ia  (Al-14)  and  (Al.iS)  is  tbsn,  ia  sa  obvious 
KCtsiioa, 


mm) = expixMmH /l.  exp\Xtim!^] 

(A1.17) 

n(9))=2i>«o»i&) 

k 

(A1.18) 

Plk)  =  2  exflX^m)!”)! 2  2  sxpIJTj  «#)/4 

(A1.19) 

We  now  recognize  the  merit  of  ^^E^s  as  providing  a  systematic  method  for  dif¬ 
ferentiation  bettt-een  alternative  hypotheses  (beliefs,  physical  mechanisms)  on  the  basis 
of  incoming  ewdence.  The  shortcoming  of  is  that  much  of  our  knowledge,  under¬ 
standing  and  ‘feeling*  for  the  behariour  of  the  physical  ^'stem  that  is  not  phrased 
explidtly  in  the  form  of  constraints  receives  no  expression  in  the  predictions. 


Appendix  2.  Fracture  dhersity 
A2.1.  Causes  of  fraaure 

Mechanical  fracture  is  frequently  the  physical  process  that  precedes  and  causes 
ftagmentation,  though  such  is  not  the  case  in  droplet  formation  (section  6.5), 
photohagmentation,  nuclear  events  (section  6.6),  desorption,  etc. 

At  the  root  cf  the  fracture  of  solids  is  the  breakage  of  atomic  or  molecular  bonds. 
Simulations  of  the  fracture  process  on  this  fundamental  level  have  been  numerous 
(Esterlihg  1979,  Latanision  and  Pickens  1983,  Paskin  et  al  1983,  Thomson  1986, 
Thomson  and  Fuller  1982).  There  are  three  considerations  that  render  this  approach 
of  limited  scope  or  excessive  difficulty.  (‘Clearly,  for  the  near  term,  there  is  not  likely 
to  be  a  first-principles  calculation  of  a  crack  tip’  (Thomson  1986}.) 

(i)  Subsequent  to  fracture  there  is  a  large-scale  relaxation,  i.e.  shift  of  particle 
positions,  in  the  solid  surrounding  the  created  discontinuity  (or  bond  breakage)  (Mai 
and  Lawn  1987).  Thus  a  presumed  single-particle  scission  is  in  truth  a  many-particle 
effect.  Owing  to  this  and  other  causes,  nominal  atomic  surface  energy  densities  (the 
energy  input  needed  to  cut  atomic  bonds  over  a  unit  surface)  exceed  observed  surface 
energy  densities  (the  quantity  that  enters  Griffith’s  expression  for  the  fracture  strength 
(Griffith  1920))  by  two  or  three  orders  of  magnitudes  (Sih  1983). 

(ii)  The  region  in  the  vicinity  of  the  crack  tip  shov/s  (even  in  brittle  solids)  plastic 
behaviour.  Such  behaviour  is  difficult  to  model  with  elastic  strings  (obeying  Hooke's 
law,  at  least  up  to  a  range,  as  is  frequent  practice  in  atomistic  simulations)  or  to 
formulate  fracture  in  the  context  of  ordinary  thermodynamics. 

(iii)  The  applied  stress  that  brings  about  breakage  does  so  not  by  initiating  and 
opening  up  a  single  void,  but  by  giving  birth  to  a  number  of  microcracks  (or  microvoids, 
in  ductile  materials),  which  then  interact  and  ultimately  merge  to  propagate  the 
macrocrack  (Curran  et  al  1987).  This  viewpoint  is  becoming  increasingly  more  firmly 
established  as  better  observational  techniques  on  the  microevent  length  scale  (<1  ^m) 
and  fast  timescales  (<10-'^  s)  are  available  (Bowling  et  al  1987,  Faber  et  al  1988).  In 
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fast  processes  (which  cause  mcst  of  the  fiagmentatioiis  treated  in  this  review)  the 
proliferation  of  microcracks  is  particularly  dedsive.  This  can  be  understood  by  the 
following  argument- 

Two  events  taking  place  in  the  solid  at  a  separation  of  develop  independently 
of  each  other,  provided  the  ratio  propagation  velocity/A/  is  smaller  than  the  relative 
strmn  rate  eje  due  to  an  externally  applied  force.  "With  propagation  or  sound  velocities 
of  the  order  of  10^  m  s~‘,  strains  £  ~  10"*  and  strain  rates  i  >  10^  s"*  (appropriate  to 
explosive  breakage),  two  cracks  at  distances  larger  than 

A/=10"’m  (A2.1) 

prefer  to  develop  independently,  rather  than  let  stress  release  operate  by  the  growth 
of  one  crack  (presumably  the  larger  one).  Likely,  the  true  distances  are  smaller  than 
that  in  (A2.1). 


A2.2.  Fracture  in  real-life  situations 

Modalities  of  fracture  in  real-life  situations  are  so  multifarious  that  either  an  all- 
embracing  treatment  or  one  that  is  based  on  a  single  viewpoint  is  likely  to  court  failure. 
Therefore,  in  the  following  brief  essay  we  emphasize  just  the  different  varieties  of 
modes  and  forms  that  fracture  can  take,  as  seen  by  several  leading  workers  in  the 
field-  The  putative  moral  from  this  exercise  is  that  in  view  of  the  wide  spectra  of 
manifestations  one  is  inevitably  thrown  back  to  a  stochastic  description,  like  that 
provided  by  mem. 

The  three  modes  of  external  stress  with  respect  to  the  fracture  plane,  shown  in 
figure  15  and  traditionally  known  as  modes  I,  II  and  III,  are  perhaps  the  least  relevant 
subdivision  of  fracture-causing  processes,  since  in  most  of  the  fragmentation  situations 
discussed  in  this  review  the  imposed  stress  appears  in  mixed  modes.  This  claim  is 
supported  by  recent  observation  of  crack  types  in  indented  ceramics  where  also  the 
variety  of  observed  cracking  sequences  has  defied  systematization  in  terms  of  the 
nature  of  the  applied  load  (Cook  and  Pharr  1990).  In  talking  about  the  initiation  of 
microdefects,  one  should  perhaps  temper  the  claim  for  an  excessive  variety  in  the 
following  sense.  Theories  of  microvoid  nucleation,  pioneered  by  Raj  and  Ashby  (1975) 
in  the  context  of  polycrystalline  metals,  apparently  lead  to  the  same,  exponential 
microflaw  size  distribution  (as  a  function  of  the  linear  flaw  dimension),  no  matter  what 
are  the  details  of  the  flaw-producing  mechanism  (Curran  et  al  1987). 

An  important  differentiation  is  between  ductile  and  brittle  fracture,  the  former 
being  characteristic  of  most  metals,  metallic  alloys  and  also  of  high  temperatures  and 
low  strain  rates,  whereas  the  latter  operates  in  rocks,  ceramics,  etc,  at  low  tem¬ 
peratures  and  in  fast  processes  (ubiquitously).  Typical  flaws  are,  respectively,  voids 
and  cracks.  A  seminal  contribution  to  fracture  processes  in  a  metal  ring  was  made  by 
Mott  (1947),  who  used  momentum  conservation  to  calculate  the  time  that  stress 
unloading  requires  to  propagate  away  from  the  point  of  fracture  until  completion  of 
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fdlure.  For  elaboration  of  this  theoiy  we  refer  to  Kipp  and  Grady  (1985).  For  brittle 
Sacture  the  corresponding  advances  are  due  to  Griffith  (1920)  and  Irwin  (1957),  who 
formulated  an  energy-balance-based  criterion  for  the  run-away  motion  of  cracks.  This 
will  happen  when  the  release  of  stress  energy  around  a  growing  crack  (regarded  as  a 
free  surface)  exceeds  the  ener^  expenditure  involved  in  the  crack  extension. 

Using  considerations  similar  to  those  given  in  the  previous  section  regarding  the 
rate  dependence  of  defect  activation,  Grady  and  Kipp  (1989)  distinguish^  between 
fracture  processes  that  are  dominated  by  the  input  stress  energy  and  those  by  the 
defect  population  density.  In  the  former  the  rate  of  fracturing  depends  on  the  loading 
intensity;  in  the  latter  strengthening  the  load  does  not  enhance  fracture,  since  the 
number  of  microdefects  is  insufficient  to  respond  to  the  added  load  and  the  energy 
gets  diverted  into  other  channels  (like  mass  motion  or  heat).  As  a  rule,  high  strain 
rates  (i)  ensure  the  creation  of  sufficiently  dense  defect  population  so  that  fracturing 
rate  becomes  energy-dominated,  while  low  e  creates  defect  bottlenecks. 

Experimentally  observed  crack-tip  shape^fall  into  three  broad  categories  that  link 
the  fracture  process  both  to  the  brittle-ductile  material  characterization  and  to  the 
emission  of  dislocations  from  and  around  the  growing  crack  (Thomson  1986).  Sharp 
knife-edged  cracks  cleave,  have  brittle  characteristics  and  displace,  but  do  not  gener¬ 
ate,  dislocations  as  they  advance;  wedge-shaped  cracks  emit  new  distinctly  oriented 
dislocations  that  advance  to  a  free  surface  as  the  crack  grows;  finally,  blunt  or  rounded 
crack  tips  emit  a  mixture  of  dislocations,  typical  of  plastic  behaviour.  The  implication 
of  dislocations  in  crack  motion  is  well  recognized  (it  is  akin  to  the  molecular  basis  of 
gas  thermodynamics),  and  the  dislocation-based  approach  has  had  its  successes  in 
relating  metal  hardness  to  granularity  (Hirth  1972).  However,  difficulties  of  treatment 
have  so  far  prevented  a  broadly  significant  contribution  to  crack  growth  and  fracture. 
Perhaps  recent  advances  in  lattice  gauge  theories  might  point  to  new  directions,  though 
at  this  stage  this  clearly  seems  a  remote  subject  on  which  to  pin  our  expectations 
(Kleiriert  1989). 
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Abstract — A  method  for  calculating  the  dependence  of  the  Gruneisen  parameters  of  a  cracked  solid  on 
crack  density  is  presented.  The  method  utilizes  the  elTective  elastic  moduli,  for  which  two  alternative 
theories  are  employed,  the  self-consutent  scheme  and  the  difTerential  scheme.  Numerical  calculations 
indicate  that  the  Grundsen  parameters  decrease  with  increasing  crack  density,  the  rate  of  decrease  being 
higher  for  the  self-consistent  scheme. 

Keywords:  Gruneisen  parameter,  cracks,  Debye  model,  effective  elastic  constants. 


1.  INTRODUCTION 

The  purpose  of  the  present  work  is  to  calculate  the 
Gruneisen  parameter  of  an  elastic  body  containing 
many  disordered  micro-cracks.  Since  cracked  solids 
-  may  be  ^ewed  as  a  particular  subset  of  composite 
materials,  we  first  recall  what  is  known  about  the 
Gruneisen  parameter  of  the  latter.  Budiansky  [1]  has 
presented  a  theory 'for  evaluating  the  Gruneisen 
parameter  y  of  composites  consisting  of  a  random 
mixture  of  a  number  of  constituents,  mth  roughly 
sphencal  particles.  His  method  was  a  natural  exten¬ 
sion  of  previous  self-consistent  calculations  of  the 
elastic  constants  of  a  composite  medium  [2].  How- 
_  ever,  his  formula  for  y  contains  parameters  which 
are  not  easily  related  to  experimentally  measurable 
quantities. 

Many  methods  for  calculating  the  effective  elastic 
moduli  of  cracked  solids  have  been  presented  in  the 
literature  (3-9|,  but  none  of  these  seems  to  have  been 
extended  for  the  purpose  of  evaluating  •/.  We  present 
here  a  simple  scheme  for  extending  any  given  theory 
for  the  effective  elastic  moduli,  so  as  to  enable  the 
calculation  of  the  dependence  of  y  on  crack  density. 
The  only  experimental  parameters  that  are  required 
are  the  elastic  moduli  of  the  uncracked  solid  and  their 
pressure  derivatives.  We  apply  this  procedure  to  two 
effective  elastic  moduli  theories,  the  self-consistent 
scheme  (SCS)  (4]  and  the  differential  scheme  (DS) 
[7, 9],  and  present  examples  of  calculations  of  y  as  a 
function  of  crack  density. 


2.  EFFECTIVE  ELASTIC  MODULI  AND  THEIR 
PRESSURE  DERIVATIVES 

We  consider  an  elastic  and  isotropic  solid  which'is 
permeated  by  many  flat  cracks.  The  distributions  of 
the  sizes,  locations  and  orientations  of  the  cracks  are 
assumed  to  be  random,  so  that  the  cracked  body  is 
macroscopically  isotropic  and  homogeneous.  The 


proposed  method  for  evaluating  y  utilizes  the  effective 
elastic  moduli  and  their  pressure  derivarives.  From 
the  many  available  methods  for  calculating  the  effec¬ 
tive  moduli  we  employ  here  two  representative  ones, 
one  of  the  SCS  type  and  one  of  the  DS  type.  While 
in  both  theories  the  effective  elastic  moduli  decrease 
with  increasing  crack  density,  the  most  obvious 
difference  between  the  two  schemes  is  that  in  SCS  the 
effective  moduli  vanish  at  a  finite  value  of  the  crack 
density  [4],  whereas  in  the  DS  they  vanish  only 
asymptotically  with  increasing  crack  density  [9].  It  is 
not  our  purpose  here  to  discuss  the  merits  and 
disadvantages  of  the  two  methods,  but  simply  to 
point  out  that  each  of  them  can  be  employed  as  a 
basis  for  the  evaluation  of  the  dependence  of  y  on 
crack  density. 

(a)  Self-consistent  scheme 
We  employ  the  well-known  version  of  Budiansky 
and  O’Connell  [4].  For  a  solid  containing  flat  circular 
cracks  the  effective  bulk  modulus  R  is  related  to  the 
bulk  modulus  K  of  the  undamaged  solid  by 


Here  v  is  the  effective  Poisson  ratio  and  e  is  the  crack 
density  parameter.  For  circular  cracks  €=tN<o’>, 
where  N  is  the  number  of  cracks  per  unit  volume,  a 
is  the  crack  radius  and  the  angular  brackets  denote  an 
average.  The  effective  shear  modulus  G  is  related  to 
Gby 


g  32(l-v-)(S-v-) 

45  (2-v) 

and  the  dependence  of  v  on  e  is  given  by 
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No  presiure  denvative  of  <  ippun  in  these  equitioiu 
because  3<13P  «•  0.  The  long  dimensions  of  a  vruck 
are  unail'ected  by  hydrostatic  preuure  applied  to  the 
eiirrouadini  medium  (tO]  and  (h«  crack  density  pur- 
amttcr  does  not  depend  on  the  width  of  the  small 
opening  between  opposite  crack  faces.  Equations 
enable  us  to  calculate  the  effective  elastic 
moduli  and  their  pressure  derivatives,  using  the  corrs- 
spoQdlng  parameters  of  the  uncncktd  solid. 

(A)  Difftrtnlial  scheme 
We  employ  the  D,S  retulrt  for  »  body  containing 
pennyshaped  cracks  (7, 9].  The  eflective  moduli  are 
even  by 

•ts  ,  /l+7^  3.  f 

(9) 
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Diflerenlialing  (7)  (9),  again  using  SkIOP  m  0,  we 
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Equations  7-13  ng.tin  provide  u  whetne  for  calcu* 
lating  the  efTcctive  elastic  moduli  and  their  pressure 
darivatives. 

i.  UKONHISEN  PARAMn'ERS 

We  first  recall  some  properiie.t  of  the  Grimeisen 
parameter  of  perfect,  uncracked  solids.  The 
Griineiien  parameter  v.  which  appears  in  the 
Mi*-Gfun«!sen  equation  of  state,  can  be  tApressed  in 
the  fomt  (1 1). 

iK 

f-TTe’  (•■*) 

where  a  is  the  volum*  coetticient  of  thermal  expan¬ 
sion,  C,  is  the  specific  hcot  at  constant  voluiiie  and  p 
Is  the  density.  Equation  (U)  is  convenient  for  avalu* 
aiing  y  from  parameicri  for  which  txpccimenial  data 
are  usually  available.  For  our  purpose  of  calculating 
the  dependence  ofy  on  oruck  density  we  find  it  inocc 
expedient  to  start  from  the  more  baste  relation  [1 1] 


which  cxpiojcs  y  in  icims  of  the  microscopic  y,  of  the 
lattice  vibraiinnxl  mruirs.  The  mode  Grunoisen  par¬ 
ameters  ate  defined  by 

i  In  V, 

where  v,-  is  the  frequency  and  c,  the  Einstein  heat 
capacity  of  the  i-ih  normal  mode.  The  summaiiona 
in  (15)  arc  taken  over  all  the  iiuiiiiai  modes  of  the 
lattice.  Calculations  of  v,  for  all  acoustic  and  optical 
modes  have  been  p'crfocmed  for  some  sulids  using 
detailed  lattice  dynamical  models  {12-14],  For 
cracked  solids  such  an  approach  will  be  iiiipmuilc- 
able  and  we  therefore  resort  to  a  Debye  type  Approxi¬ 
mation,  in  which  the  summations  are  taken  only  over 
the  three  acoustic  hnmehes.  which  are  assumiid  to  be 
non-difpersivc.  For  an  acousu'c  mode  chatuuierizcd 
by  an  clasilc  constant  Cl  the  mode  GruneiMn  par¬ 
ameter  is  given  by  (he  Smith  formula  [20] 

1  K  CC, 

V,  m  — •  +  •— , 

'*  fi  2C,  (5/ 


(17) 
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For  an  isotropic  solid  there  are  only  two  sound 
velocities,  the  shear  velocity  v„  given  by  pv;  =  G, 
and  the  compressional  velocity  given  by 
pvl  =  K  +  4GI3.  The  corresponding  Gruneisen  par¬ 
ameters  are  given  by: 


At  very  low  temperatures  only  the  non-dispersive 
acoustic  modes  are  excited  and  the  specific  heat 
weighting  in  eqn  (15)  may  be  replaces  by  a  v~^ 
weighting,  where  v  is  the  wave  velocity.  The  low 
teiuperature  limit,  •/(,■  of  the  Gruneisen  parameter  is 
therefore  ^ven  by 

•/L  =  7X-(“^yr  +  2y,),  (20) 

where  a  =  vjv^.  For  the  high  temperature  Gruneisen 
parameter  Vh  ihe  Debye  model  yields  the  simple 
average 

7H  =  i(yp  +  2y.).  (21) 

While  expression  (20)  for  '/l 's  exact,  eqn  (21)  for  Vh 
is  an  approximation.  Nevertheless,  for  a  number  of 
materials  '/h  lies  close  to  the  room  temperature  value 
of  y,  as  determined  experimentally  from  (14). 

So  far  we  have  discussed  the  Gruneisen  parameters 
of  the  perfect  solid  (e  ^  0).  The  dependence  of  the 
Gruneisen  parameters  on  crack  density  is  obtained  by 
replacing  k,G,  dK/dP  and  dG/dP  in  all  the  formulae 
of  this  section  by  K,  C,  dR/dP  and  dGjdP,  respect¬ 
ively.  The  latter  quantities  arc  calculated  by  the 
methods  de.scribcd  in  section  2. 


4.  NUMERICAL  RESULTS  AND  DISCUSSION 

We  present  the  results  of  calculations  of  the  depen¬ 
dence  of  the  Gruneisen  parameters  on  crack  density 
for  the  minerals  periclasc  and  lime.  They  were  chosen 
as  examples  because  for  £  =  0  their  calculated  high 
temperature  Gruneisen  parameters  Vh  agree  well  with 
the  experimental  room  temperature  values  obtained 
from  (14)  (21),  Thus,  we  may  expect  •/«  to  provide  a 
good  estimate  of  •/  also  for  the  cracked  solid.  For  the 
uncracked  solid  elastic  moduli  and  their  pressure 
derivatives  the  data  compiled  by  Anderson  er  al.  (21) 
have  been  employed.  For  each  material  both  SCS 
and  DS  calculations  were  performed.  The  calculated 


Fig.  I.  Calculated  dependence  of  the  effective  bulk  modulus 
on  crack  density  for:  (a)  periclase;  (b)  lime.  Full  curves — 
SCS;  dashed  curves — DS. 


dependence  of  the  effective  moduli  R  and  G  on  crack 
density  is  shown  in  Figs  1  and  2,  respectively.  The 
Gruneisen  parameters  are  shown  in  Figs  3-6.  We  find 
that  the  Gruneisen  parameters  decrease  with  increas¬ 
ing  crack  density.  The  rate  of  decrease  is  higher  in  the 
SCS  than  in  the  DS.  In  this  respect  the  qualitative 
benaviour  is  similar  to  that  of  the  effective  elastic 
moduli  (Figs  I  and  2),  which  also  decrease  faster  in 
the  SCS. 

Acknowledgements— work  was  sponsored  in  part  by 
U.S.  Air  Force  Systems  Command,  US.AF  under  Grant  No. 
89-0374. 


Fig.  2.  Calculated  dependence  of  the  effeettve  shear  modulus 
on  crack  density  for;  (a)  periclase;  (b)  lime.  Full  curves — 
SCS;  dashed  curves — DS. 


Fig.  4.  Dependence  of  Gruneisen  parameters  of  penclase  on  Fig.  6.  Dependence  of  Gruneisen  parameters  of  lime  on 
crack  density,  as  calculated  from  the  DS.  crack  density,  as  calculated  from  the  DS. 
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R.  Englmaa: 

’’Statistical  Models  for  the  Fracture  of  Disordered  Media”  Edited  by  H.J. 
Herrmann  and  S.  Roux  (North  Holland,  Amsterdam,  1990) 

A  Review  for  J.  Statistical  Physics. 

The  book  divides  easily  into  three  parts:  An  introductory  one,  an 
experimental  part  (or,  more  accurately,  a  material  science  part)  and, 
finally,  the  detailed  description  of  several  theoretical  models. 

Most  of  the  introductory  chapters  are  written  by  the  editors  and  are  a  joy 
to  read:  especially  for  being  shown  how  engineering  concepts  are  moulded 
into  physical  thinking  or  language  (like  thermodynamic  energy,  duals, 
localization,  renormalization,  scaling,  etc.) 

The  experimental  sections  focus  on  the  dependencies  of  material  properties 
on  parameters  like  temperature,  composition,  morphology,  manufacture, 
treatment,  changes  in  microstructure,  process  rates. 

When  one  reaches  the  theoretical  papers  one  is  confronted  with  the 
realisation  that  these  hardiy  ever  address  the  real  problems  and  the 
questions  that  occupy  practitioners.  (In  essence,  these  questions  add  up  to 
one:  how  to  produce  materials  that  have  improved  material  properties.  In  the 
ceramic  -  fish  diagram  on  page  44  each  scale  stands  for  a  real-life  problem 
that  awaits  solution.  Needless  to  say  that  scaling-laws  do  not  provide  this. 
An  exception  is  the  article  on  fragmentation,  which  is  geared  to  coramunition 
processes,  though  only  in  a  qualitative  fashion. However ,  in  another  article 
(on  randomness)  one  experimental  figure,  about  twenty  years  old,  is  all  that 
is  presented'  in  ’’support"  of  25  pages  of  results  from  models. 

It  is  clear  that  fracture  is  a  difficult  subject  whose  complexities  have 
defied  70  years  of  research.  The  approach  in  this  book  (and  the  uniformity 
of  the  articles  bemuses  the  reader)  breaks  away  from  the  traditional,  for 
reasons  that  (I  would  hazard)  are  rooted  in  the  shine  and  lure  of  moderm 
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has  so  few  experimental  facts  for  its  s.upport,  a  more  balanced  choice  of  the 
articles  would  have  been  justified  in  a  book  of  this  title.  Both 
’’Statistical  models"  and  "Disorder”  are  also  present  in  the  approaches  of 
Atkins,  Evans,  Kachanov,  Lawn,  Mai,  Rice,  where  disorder  enters  in  the  form 
of  microcracks,  voids,  stress  distributions  either  before  or  after  the 
application  of  external  stresses. 

In  many  of  the  simulations  described  in  the  book  there  are  springs  or 
beams  (that  oppose  fracture  or  bending  on  a  basic  level),  and  it  is  not 
clear  how  these  relate  to  the  primary  bonds  that  characterize  the  material 
The  question  is  not  pedantic,  since  it  is  these  primary  bonds  that  the 
material  scientist  aims  at  changing  (by  introduction  of  new 
constituents,  additives  or  reinforcements)  to  meet  his  ends.  Similarly,  it 
is  not  evident  how  the  primary  length  unit  in  the  simulation  is  chosen  (in 
cases  other  than  granular  solids,  in  which  the  grain  may,  but  not  necessarily 
does,  define  a  fundamental  length  scale).  One  supposes  that  an  averaging  of 
some  material  property  over  some  domain  in  the  medium  is  implied  and  that  the 
size  of  this  domain  defines  the  basic  length  scale.  But  this  is  not  discussed 
and  one  wonders  in  view  of  instability  of  some  cracks  whether  the  length  unit 
ought  not  to  depend  on  stress  or  on  some  other  parameter. 

A  quotation  from  page  150  in  the  book  well  summarizes  the  situation: 

"The  authors  have’’"presented  their  belief  of  what  the  Truth  is,  but  the 
reader  should  keep  in  mind  that  there  are  other  competing  Truths  around.  ’’ 
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PART  8 

ggggturg  in  granular  materials  with  oeoraetrlcal  disorder 

(Work  performed  by  M.  Lemanska,  Z.  Jaeger  and  R.  Englman) 

A  two  dimensional  model  of  solids  consisting  of  polygonal  grains  has 
been  constructed  by  Voronoi-tessellation.  This  means  that  (up  to  2000) 
points  were  chosen  randomly  in  a  plane  and,  by  joining  up  the  perpendicular 
bisectors  that  link  neighboring  points,  the  plane  was  divided  into  polygons 
of  different  sizes  and  shapes  (Figure  1) .  The  average  number  of  sides  in 
the  polygons  is  6.  For  a  theory  of  side  distribution  in  2D  and  in  3D  see 
Ref.  1.  The  present,  Voronoi-tessellation  is  a  reasonable  approximation  to 
disordered  multiphase  structures  containing  irregular  grains.  More 
precisely,  it  is  deemed  superior  to  alternative  models  in  which  disorder 
enters  through  randomness  of  bond  strengths  in  a  geometrically  regular  (e.g. 
square  or  triangular)  lattice. 

Fracture  is  simulated  by  applying  either  a  fixed  stress  or  a  fixed 
elongation  (constant  strain)  at  two  parallel  edges  and  calculating  the 
elongation  in  each  bond  for  a  harmonic,  central  force  model  of  springs. 

One  or  more  bonds  break  when  their  calculated  extensions  or  relative 
extension,  (defined  as  the  ratio  of  the  extension  to  the  original  bond 
length)  exceed  a  certain  value.  Then  the  stress  or  strain  is  reapplied  or 
is  raised  until  newer  bonds  breaks  (and  in  fact  they  do  so  with  greater 
facility  since  the  broken  bonds  no  longer  support  the  stress) .  The 
procedure  is  repeated  until  one  gets  a  continuous  domain  of  broken  bonds 
cutting  across  the  sample,  which  is  then  considered  fractured. 

Several  parameters  of  the  fracture  process  are  computed,  such  as: 

(a)  the  least  stress  or  strain  v/hich  causes  fracture  (the  fracture  stress 
or  strain) , 

(b)  The  number  of  bonds  that  are  broken  at  each  stage  of  the  fracture 
process, 

(c)  the  location  of  the  broken  bonds  (Fig.  2) 

(d)  the  distribution  of  broken  bonds  as  function  of  their  length 
(Fig.  3). 
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(e)  tise  ursjer  of  bssasJs  broiea  af  zc^ttsre  ss  ftnicfica  of  sasrolc  S2.ze 
CFig.  4)  asd  of  ctfeer  oarsaebers- 

Sie  work  is  still  ia  progress,  witb  a  view  of  gaining  insigSst  into  tJtc 
conaectioa  between  sicrepreperties  and  breakdown  strengtia.  Still,  cae  can 
already  derive  Fig.  4-  an  iaqportant  connecticn  between  tfce  present 

Vorcr-oi  anH  »-;w  mnericaliy  isore  extensively  studied  regular 

lattice  zsodels  with  streagtb  disorder: 

ia  tee  latter  at  breakdown,  tiie  nirxyer  of  fractured  bends  n^  is  found 
to  be  related  to  the  linear  size  L  of  the  ^stea  by: 

(Ref.  (3)) 

The  slope  of  1.85  ia  Fig.  4  is  quite  close  to  this,  shewing  the 
similarity  in  the  breakdown  behavior  of  solids  that  are  elastically  and 
geometrically  disordered. 
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Figure  C«pti<m» 


Fig.  1.  Voronoi-tessellatioa 

A  relatively  saall  aua&er  of  points  (70,  aisEbered  1-70)  are  placed 
randoaly  in  the  plane  and  the  area  sub-divided  into  polygons  {as 
described  in  the  tert) . 

Fig.  2.  Fracture  in  a  disordered  structure. 

In  the  last  of  three  steps  (achieved  by  successively  reapplying 
the  strain)  the  structure  is  ruptured  by  intragranular  cracking. 
Full  lines  are  normal  to  intact  intragranular  surface.  Broken 
lines  represent  cracked  surfaces  normal  to  the  lines. 

Fig.  3.  Histogram  of  broken  bonds  as  function  of  the  original  bond  length 
at  fracture.  The  data  represent  averaging  over  several  runs.  The 
criterion  for  bond  rupture  is  that  the  relative  elongation  6  of 
bonds  exceeds  0.2.  The  distribution  of  original  lengths  is  also 
shown  by  broken  lines  and  the  percentage  of  bonds  broken  by  dotted 
lines.  In  our  model  short  bonds  have  a  preference  for  being 
broken . 


Fig.  4.  Log-log  plot  of  the  number  of  broken  bonds  (at  rupture  of  the 
sample)  versus  the  square  root  of  the  number  of  Voronoi  polygons 
(essentially,  the  linear  size  of  the  system) 

The  slope  of  the  line  is  about  1.85. 
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PART  9 


Experiawntal  Validation 

Extensive  tests  of  fracture  and  fragmentation  behaviors  of  reinforced 
concrete  slabs  have  been  performed  by  Vargas,  Hokanson  and  Rindner(l). 
Since  their  test  results  are  readily  interpretable  in  terms  of  our 
percolation  approach  (especially  through  their  employment  of  varicible  length 
scales  for  reinforced  bar  spacing  and  slab  thickness) ,  we  describe  the 
experimental  arrangement  at  some  length. 

Concrete  walls  with  dimensions  (in  cm  length  units)  of  about  50  x  50  x 
(5-8)  were  subjected  to  detonations  at  close  quarters  from  several  types  of 
charges.  The  walls  were  reinforced  by  tightly  locking  horizontal  and 
vertical  wires  in  the  concrete  matrix,  with  spacing  that  varied  from  2.5  to 
5  cm.  The  debris  was  collected  and  its  distribution  determined  as  function 
of  mass,  velocity,  range  and,  in  a  less  precise  manner,  of  shape  and 
provenance.  The  present  section  treats  only  the  mass  distributions,  as 
being  most  closely  related  to  the  theories  laid  out  in  the  preceding  parts; 
the  analysis  of  the  other  distribution  requires  more  detailed  considerations 
of  the  dynamics  of  the  fragmentation  process,  which  has  not  yet  been  given. 

In  summary  it  can  be  -tated  that,  though  the  experimental  data  are 
subject  to  a  great  deal  of  scatter,  there  are  strong  indications  of  a 
percolation-like  process  occurring  causing  the  fragmentation. 

Following  the  authors  of  the  report  ,  we  quantify  the  external  forces 
that  bring  about  the  fragment  by  the  specific  impulse  factor  IF  which  is 
defined  as  "the  total  impulse  (imparted  by  the  exploding  charge  to  the  slab] 
divided  by  the  square  root  of  the  effective  [slab]  thickness".  Without 
undertaking  a  deep  analysis  of  IF,  we  can  appreciate  that  it  affects  both 
the  strength  of  the  stresses  in  the  fragmenting  slab  and  the  strain  rate  in 
it . 

The  geometrical  parameters  are  the  reinforced-bar  ("rebar"  spacing  (Rg) 
and  the  slab  thickness  T.  The  effect  of  the  former  on  fragmentation  is  much 
more  marked  than  that  of  the  latter  (as  will  be  shown  shortly) .  This  leads 
us  to  the  supposition  (which  actually  has  some  far-reaching  practical 
implications) ,  that  the  rebars  divide  up  the  slab  into  smaller  panels  of 
size  Rg  X  Rg.  The  different  panels  fragment  independently  of  each  other;  in 
other  words;  the  effective  lateral  dimension  of  the  solid  is  Rg  (rather  than 
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the  slab  size)  .  Furthermore,  there  are  indications  from  the  data  that  the 
slab  (or  the  panels)  behaves  as  a  two-dimensional  solid  (more  precisely,  the 
effective  dimensionality  of  the  slab  is  nearer  to  2D  than  to  3D,  as  will  be 
seen  later . )  Since  the  rebars  lie  in  two  horizontal  planes  that  are  0.6  cm 
from  the  exterior  planes  and  4-6  cm  from  each  other,  it  is  not  clear  a 
priori  why  the  rebars  should  have  such  a  strong  effect  on  the  fragmentation 
process.  If  a  future  investigation  should  find  the  answer  to  this  question, 
this  could  improve  the  use  of  rebars  in  concrete  (and  also  of  whisJcers  and 
fibers  in  ceramics,  etc.). 

In  a  more  recent  worlc  on  explosion  effects  on  aircraft  shelters  (ACS) , 
fragment  distributions  (adding  up  to  40  tons)  were  obtained (for  a 
graphical  representation  of  the  results  see  (3) ) .  Though  the  test  results 
available  in  (3)  were  fewer  than  in  the  smaller  structures  (slabs)  of  (1) , 
we  have  also  included  these  in  our  percolation  theory  analysis. 

1.  Cummulatlve  mass  distributions 

Percolation  theory predicts  the  following  density  of  fragments 
(no.  per  unit  mass)  as  function  of  the  fragment  mass. 

^  eC 

where  x  =  2.05  in  two  dimensions  (2D) 

=2.2  in  three  dimensions  (3D) 

The  cumulative  distribution  will  threfore  have  exponents  -  (x  -  1)  and 
will  appear  on  a  log  -  log  m  plot  as  a  straight  line.  About  14  of  the  39 
test  results  in  the  slab  test^^^  have  a  linear  section  in  the  higher  mass 
range;  in  25  no  linear  section  was  discerned.  Examples  of  the  linear 
section  are  seen  in  Fig.  1.  The  slopes  [identified  with  -  (x  -  1)  ]  differ 
from  case  to  case.  The  cummulative  distribution  of  the  14  slopes  is  shown 
in  Fig.  2.  Though  a  significant  scatter  is  evident  among  the  test  events, 
the  slopes  fall  quite  close  to  the  percolation  theory  values :  1.05  (two 
dim.)  and  1.2  (three  dim.),  showing  better  agreement  with  the  former.  We 
have  also  included  the  two  slopes  from  ACS  results ,  which  also  lie  near 
the  percolation  -  theory  slopes. 
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Let  it  be  remarked  that  the  linearity  in  the  long-log  plot  by  itself  is 
a  strong  vindication  of  the  percolation  model,  since  most  of  the  other 
distributions  are  exponential .  The  curving  (non-linearity)  of  the 

distributions  in  the  smaller  mass  range  is  an  indication  of  some  intrinsic 
length  scale,  like  some  structure  or  granularity  of  the  material. 

II.  Impulse  dependence  of  fragments 

The  authors  of  the  slab  tests counsel  a  more  detailed  statistical 
investigation.  This  was  undertaken  by  us  for  both  the  observed  total 

fragment  number  N  and  mean  fragment  mass  m  collected  as  functions  of  the 

impulse  factor  IF.  Following  Ref.  (1)  we  have  made  a  fit  with  linear 

regression  of  the  form: 

N  -  A(IF)  +  B 

m  =  C(IF)  +  D 

separately  for  four  combinations  of  the  geometrical  parameters  slab 
thickness  T  and  rebar  spacing  Rg.  The  results  are  shown  in  Table  1.  The 
regressions  were  examined  for  statistical  significance  by  Student's 
t-test^®^  .  The  obtained  t-values  (shown  in  columns  7  &  12  of  the  table)  are 
much  higher  (better)  than  those  (2-4)  needed  for  assuring  significance.  At 
this  stage  the  regression  analysis  does  not  warrant  the  drawing  of  physical 
conclusions  -  except  from  the  mean  mass  m  (last  column  in  the  table) . 

Table  1.  Statistical  Analysis  of  Impulse  Dependence  for  Fragment  Number  and 

Size  (in  cm) .  (Mass  in  grams) . 


Geometry 

T  Rg 

#  Event 

n 

Regression  coefficients  and  significance 

A 

B 

a 

i 

N 

C 

D 

B 

B 

m  1 

5 

2.5 

9 

12.5 

123 

.09 

25.5 

167 

-5 

31 

.13 

18 

13 

5 

5 

8 

36 

36 

.05 

47 

180 

-1.8 

45 

.17 

13 

3 

7.5 

2.5 

9 

2.7 

88 

.14 

16 

100 

-1.4 

21 

.12 

14 

15 

7.5 

5 

9 

35 

-6 

.09 

25 

146 

-19 

-37 

.08 

20 

45 
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III.  Finite  size  scaling  of  inean  niasg. 

The  entries  for  m  show  that  the  mean  mass  scales  strongly  with  rebar 
spacing 

m  ^  R3l-« 

and  weakly  with  slab  thickness 

m  tC  tO-^2 
(Fig.  3) 

Regarding  the  thickness  as  constant,  the  upper  relation  can  be 
rewritten  in  terms  of  the  volume  V  of  the  panel,  since  V  = 

m  «C  V  (T  constant) 

Now  the  mean  mass  of  percolating  clusters  at  the  critical  crack  density 
(p  =  p^)  scales  with  the  finite  size  of  the  sample  as 

m  cC  V(2°f  " 

where  Df  is  the  fractal  dimension  of  the  crack  cluster  and  d  is  the 
Eucledian  dimension.  (Table  II)  . 

Table  II.  Percolation  Exponents 


I - 1 - 1 

I  Dimension  d  1  2  ^  1 

I  I  1 

I  Fractal  dim.  |  1.9  2.5  | 

I  1  I 

I  2Df/d  -  1  j  .9  .66  I 

I  i  I 

I  Experiment  I  .815  I 

I _ I _ I 
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From  the  table  one  notes  that  the  experimental  exponent  lies  between 
the  percolation  exponents  in  two  and  three  dimension,  favoring  (as  in  the 
earlier  instance)  a  two  dimensional  crack  framework.  The  observed  exponents 
differ  significantly  from 

T  -  1  =  .67 

which  is  known  from  the  Gaudin  -  Schuhmann  distribution (7, 8) .  This 
characterizes  several  mineralogical  fragments  in  the  low  mass  range,  while 
our  power  law  is  for  the  heavier  masses  . 

Visual  (microphotographic),  as  well  as  quantitative,  support  for  the 
percolation  concept  is  obtained  from  creep  in  several  structural  ceramics. 
Wiederhorn  and  coworkers  have  discerned  cavity-coalescence  induced  failure 
in  some  of  their  materials  (like  Coors  SCRB210)  and  have  similarly 
interpreted  results  by  others .  The  percolation  process  observed  in  void- 
coalescence  is  highly  correlated,  anisotropic  and  has  a  position  dependence 
in  terms  of  the  distance  from  the  main  crack.  All  these  ramifications  are 
present  in  our  general  thermo-dynamic  theory  and  in  some  of  the  computations 
presented  in  this  report .  However,  for  a  quantitative  comparison  with 
available  creep  studies  and  the  derivation  of  general  laws  further  work  is 
needed. 
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Figaxe  C^ptioss 


Fig.  1.  Ccazalatlve  zsss  distribtrtion  in  coarete  slab  tests^~^.  The 
double  logarlfcinsic  plot  saovs  in  tiis  hicber  mss  range  a  linear 
portion,  indicating  an  inverse  power  law  relation. 

Fig.  2.  Distribation  of  mss  exponents.  The  cussnlative  nnnber  of  the 
observed  mss  exponents  is  shown  by  full  lines  for  the  slab 
tests  and  by  thickened  lines  for  the  aircraft  shelter  (hCS) 
tests.  Szperinental  results  show  a  scatter  around  (and  not  far) 
fron  the  percolation  theory  predictions  (broken  line  in  three 
dinensions,  chained  line  in  rwo  dimensions) . 

Fig.  3.  The  average  msses  n  against  the  slab  thickness  T  and  the  rebar 
spacing  Rg  on  a  log-log  plot. 

Full  squares  S  vs  Rg 
Open  circles  m  vs  T 

The  full  line  gives  the  exponent  1.63  of  Rg,  and  the  two  broken 
lines  both  fit  an  exponent  of  0.42  for  T. 
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Fig.  2 
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